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AN  INTRODUCTION  TO  ARRSTATS  —  A  COMPUTER 
PROGRAM  FOR  SIMULATING  THE  EFFECTS  OF  ERRORS  IN 
TIME-  AND  PHASE-STEERED  PLANAR  ARRAY  ANTENNAS 


1.  INTRODUCTION 

Phased  array  antennas  are  remarkable  for  their  suitability  to  many  applications,  which  is  partly 
because  they  steer  quickly,  allow  adaptive  processing,  conform  to  special  shapes,  and  produce  a  variety  of 
radiation  patterns.  As  low-sidelobe  radiation  patterns  were  being  developed,  it  became  apparent  that 
errors  in  the  phase  or  amplitude  of  the  element  excitations  would  limit  the  lowest  achievable  sidelobe 
level.  Such  errors  might  be  due  to  manufacturing  tolerances  on  the  physical  structure,  inaccurate  phase 
shifters,  or  non-uniform  feeds,  and  can  affect  the  average  and  peak  sidelobe  levels,  beam  width,  pointing 
accuracy,  and  directivity,  for  example.  Although  one  cannot  predict  the  performance  of  a  given  antenna 
without  measurements  of  the  phase  and  amplitude  accuracy  of  each  element,  one  can  draw  conclusions 
about  the  behavior  of  an  ensemble  of  antennas  that  are  statistically  identical  but  have  different  realizations 
of  phase  and  amplitude  errors.  The  present  work  applies  this  approach  to  a  generalized  array  antenna  that 
blends  phase-  and  time-steering  to  achieve  greater  bandwidth  at  a  reasonable  cost.1-3  The  purpose  of  this 
report  is  to  introduce  the  reader  to  a  computer  program  for  simulating  time-  and  phase-steered  planar 
array  antennas  subject  to  deterministic  and  random  excitation  errors.  It  is  intended  as  an  overview  to  the 
features  and  capabilities  of  the  program  and  a  guide  to  understanding  the  program’s  input,  processing,  and 
output. 


The  theoretical  study  of  errors  in  array  antennas  has  produced  a  large  body  of  literature. 
Introductions  to  these  works  and  key  results  and  derivations  may  be  found  in  several  books.4-6  The 
literature  on  simulations  is  more  sparse,  but  two  programs  have  been  described  recently  against  which  this 
work  may  be  contrasted.  First,  Chrisman7  describes  a  program  that  simulates  phase-steered  planar  arrays. 
It  computes  the  directivity  and  cuts  of  the  design  and  expected  radiation  patterns  from  the  error  statistics 
using  theoretical  formulas.  One  pattern  cut  intersects  the  main  beam  and  boresight,  and  the  other  is 
normal  to  it  through  the  main  beam;  from  them  the  beam  width  is  obtained  in  those  directions.  Second, 
Wright8  briefly  discusses  the  features  of  a  simulation  of  phase-steered  arrays.  Its  parameters  include  phase 
and  amplitude  errors,  number  of  quantization  bits,  and  bandwidth,  and  it  can  output  two-dimensional 
beam  patterns,  sidelobe  statistics,  and  measures  of  specialized  interest.  The  program  discussed  here, 
called  ARRSTATS,  differs  from  those  in  Refs.  7  and  8  primarily  in  that  it  can  model  arrays  with  a  hybrid 
phase-  and  time-steered  architecture,9  including  strictly  phased  arrays  and  strictly  time-steered  arrays. 
Also,  it  computes  individual  realizations  of  the  hemispherical  radiation  pattern  and  obtains  most  measures 
of  the  array  performance  directly  from  the  full  pattern  rather  than  from  formulas  or  from  pattern  cuts 
chosen  a  priori.  (Ref.  8  does  not  specify  how  its  performance  measures  are  obtained.)  For  example,  the 
beam  width  cuts  are  always  along  the  major  and  minor  axes  of  the  beam  width  ellipse,  regardless  of  its 
orientation.  Another  special  feature  is  a  method  for  locating  the  beam  peak  that  is  highly  accurate  for 
nearly  flat  phase  fronts;  it  is  the  only  measure  not  obtained  directly  from  the  radiation  pattern. 

In  more  detail,  the  modeled  array  comprises  a  rectangular  grid  of  phase-steered  elements;  these 
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are  grouped  into  time-steered  subarrays,  which  in  turn  are  grouped  into  time-steered  subapertures.  The 
phases  and  times  may  be  quantized.  A  random  phase  error  may  be  associated  with  the  elements,  and 
random  time  errors  may  be  assigned  to  the  subarrays  and  subapertures.  Also,  random  amplitude  errors 
may  be  associated  with  the  elements,  subarrays,  and  subapertures.  The  user  specifies  the  design 
frequency,  at  which  the  error-free  times  and  phases  would  correctly  steer  the  antenna,  and  the  operating 
frequency,  at  which  the  array’s  behavior  is  simulated.  The  antenna  may  be  steered  to  any  direction. 
ARRSTATS  assesses  the  array’s  performance  by  computing  and  analyzing  the  far-field  radiation  pattern  or 
an  ensemble  of  statistically  identical  patterns.  It  ignores  polarization  and  mutual  coupling  between 
elements  and  assumes  that  the  elements  radiate  uniformly  into  the  forward  hemisphere.  For  each 
computed  pattern,  it  determines  the  following: 


the  location  and  power  density  of  the  main  beam  peak 
the  error  in  the  main  beam’s  location 

the  main  beam’s  angular  limits,  major  and  minor  widths  at  half  power,  and  orientation 
the  directivity 

the  ratios  of  powers  in  the  main  beam,  sidelobes,  and  the  radiation  hemisphere 

the  powers  and  distances  of  the  sidelobes  that  are  strongest  and  nearest  to  the  main  beam 


Furthermore,  it  can  track  these  measures  as  functions  of  a  user-specified  independent  variable.  ARRSTATS 
provides  textual  output  of  the  performance  measures,  plots  of  radiation  patterns,  and  plots  of  performance . 
measures  versus  the  independent  variable.  •.  >  -T,r.  br® 

*  .  ’  ,,  ,  .  .J.-; 

ARRSTATS  is  a  script  written  for  version  5.3  of  MATLAB,  a  commercial  software  package  for. 
technical  computing.11  Some  elements  of  the  program  structure  and  some  of  the  graphics  facilities  take 
advantage  of  features  in  version  5.3,  but  much  of  the  code  is  compatible  with  earlier  versions  of 
MATLAB.  ARRSTATS  consists  of  one  text  file  and  is  executed  by  typing  its  filename  at  the  MATLAB 
command  prompt.  The  program  does  not  have  an  input  user  interface;  instead,  the  user  hardcodes  input 
into  the  script  before  execution.  These  input  points  are  tagged  with  the  word  “INPUT”  in  the  code’s 
comments. 


The  remainder  of  this  report  is  structured  as  follows.  Section  2  introduces  the  coordinate  systems 
used  for  input  and  output.  Section  3  explains  the  input  parameters  that  describe  the  array.  Section  4 
specifies  the  model  for  the  excitations,  and  Section  5  outlines  the  calculation  of  the  radiation  pattern. 
Section  6  describes  the  pattern  measures.  Section  7  discusses  looping  over  multiple  realizations  and 
parameter  values,  and  Section  8  exhibits  the  program’s  output.  Two  appendices  are  provided:  Appendix  A 
relates  the  variables  used  in  this  report  and  in  the  program,  and  Appendix  B  is  a  listing  of  the  program 
code.  Most  of  the  report  aims  to  describe  aspects  of  the  program’s  operation  but  does  not  give  the  details 
of  the  implementation  or  algorithms.  The  interested  reader  is  directed  to  the  code  listing,  specifically  the 
comments  that  introduce  each  section  of  the  program.  Throughout  this  report,  code  variables  are  printed 
in  a  monospaced  face,  and  braces  ({})  enclose  references  to  code  line  numbers  except  where  the 
context  suggests  set  notation. 


2.  COORDINATE  SYSTEMS  AND  PROJECTIONS 


ARRSTATS  internally  uses  a  three-dimensional  Cartesian  coordinate  system  to  describe  space.  The 
array  lies  in  the  x-y  plane  and  radiates  into  the  half-space  z  >  0,  as  in  Figure  1 .  (Although  we  describe  a 
transmitter  array,  the  case  for  a  receiver  array  is  identical.)  The  far-field  spatial  distribution  of  this 
radiation  —  that  is,  the  radiation  pattern  —  is  a  function  of  direction  in  the  half-space.  Two  variables 
suffice  to  specify  direction,  and  several  pairs  of  variables  are  useful  for  this  purpose.  First,  direction 
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Figure  1  —  Cartesian  and  spherical  coordinates;  6,  ty,  a,  and  e  are  shown  positive 


cosines  are  the  natural  coordinates  for  calculating  the  array  factor,  as  will  be  seen  later.  The  direction 
cosines  for  a  given  direction  are  simply  the  Cartesian  coordinates  (£  tj,  £)  of  the  corresponding  unit 
vector.  Specifying  a  direction  in  the  half-space  z  >  0  requires  only  and  tj;  £  may  be  obtained  from  the 
unit  vector  constraint  if  needed.  Second,  spherical  coordinates  are  convenient  for  constructing  flat 
projections  of  the  radiation  pattern,  as  it  may  be  regarded  as  a  function  of  location  on  a  (curved) 
hemisphere.  As  shown  in  Figure  1,  the  polar  angle  6  for  a  given  vector  is  the  angle  between  the  positive  z 
axis  and  the  vector,  while  the  azimuth  angle  y/  is  the  angle  between  the  positive  x  axis  and  the  projection 
of  the  vector  onto  the  x-y  plane.  Third,  traditional  antenna  coordinates  connect  these  simulations  to  an 
established  context.  Given  the  projection  of  a  vector  onto  the  x-z  plane,  the  azimuth  angle  a  is  the  angle 
between  the  projection  and  the  z  axis,  while  the  elevation  angle  e  is  the  angle  between  the  projection  and 
the  vector.  These  three  sets  of  coordinates  are  related  according  to 

£  =  sin#cos^  =-cosesina 
7=sin0sin^=  sine 
C  =cos0  =  cosecos  a 

cosO  =  cosecos  a  =C 
tant/A  =  -tane/sina  =tjl £ 

-tan a  ~^/C  =  tan#cos^ 
sine  =  ?/  =  sin#  sin  y/. 

ARRSTATS  also  employs  three  projections  of  the  hemisphere  onto  flat  two-dimensional  space: 
orthographic,  Lambert  azimuthal,  and  stereographic.  The  orthographic  projection  yields  a  view  of  the 
hemisphere  from  a  particular  vantage  point  without  perspective  distortion  and  is  available  for  displaying 
the  radiation  pattern.  The  remaining  projections  both  map  the  hemisphere  to  a  disk  such  that  boresight 


(1) 

(2) 

(3) 
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corresponds  to  the  center  of  the  disk  and  grazing  directions  correspond  to  the  perimeter  of  the  disk.  To 
describe  these  projections  more  specifically,  we  denote  locations  on  the  disk  using  polar  coordinates 
(radius  and  angle).  For  both  projections,  the  angular  coordinate  is  set  equal  to  the  spherical  azimuth  angle 
if/;  the  mapping  from  the  spherical  polar  angle  6  to  the  radius  r  distinguishes  the  two  projections. 

The  Lambert  projection  preserves  relative  area:  the  ratio  of  two  areas  on  the  hemisphere  equals 
that  of  the  corresponding  areas  on  the  projected  disk.12  This  property  may  be  expressed  by  equating  (to  a 
proportionality  constant)  the  spherical  and  planar  surface  areas: 

sin&dOdy/  =  CrL  drL  dy/ ,  (4) 

where  the  subscript  “L”  distinguishes  the  radius  in  the  Lambert  projection  from  that  in  the  stereographic 
projection  below.  Canceling  dy/  and  integrating  both  sides  produces  an  integration  constant,  whose  value 
and  that  of  C  are  determined  by  the  constraints  that  rL  =  0  when  6  -  0  and  rL  =  1  when  6  =  nl  2.  One  finds 

rL  =  V2  sin— ,  (5) 

L  2 

which  deviates  only  slightly  from  a  linear  relationship  for  6  e  [0,  k  /  2].  When  the  radiation  pattern  is 
plotted  with  the  Lambert  projection,  the  areas  occupied  by  structures  such  as  the  main  beam  and  sidelobes 
are  in  true  proportion  to  each  other  and  to  the  total  area. 

While  the  stereographic  projection  does  not  preserve  area,  it  is  conformal.12  The  local  scale  is 
uniform  in  any  direction;  shape  is  preserved  locally.  Great  and  small  circles  on  the  hemisphere  are 
projected  into  circles  or  straight  lines,  and  the  angle  between  two  great  circles  on  the  hemisphere  equals 
the  angle  at  the  intersection  of  their  projected  images.  To  derive  the  governing  relationship,  equate  the 
aspect  ratios  of  orthogonal  derivatives,  as 


sin  Ody/  rsdy/ 

where  the  subscript  “S”  denotes  the  stereographic  projection.  Canceling  and  integrating  as  before  yields 


r<.  =  tan—. 
s  2 


ARRSTATS  internally  uses  the  stereographic  projection  when  identifying  the  major  and  minor  axes  of  the 
beamwidth  ellipse  (see  Section  6.3)  and  makes  it  available  for  plotting  the  radiation  pattern. 

3.  ARRAY  PARAMETERS 

Several  parameters  describe  the  antenna’s  geometry  and  associated  frequencies  {29-94}.  As 
Figure  2  illustrates,  the  antenna  elements  occupy  a  regular  rectangular  grid  in  the  x-y  plane;  the  element 
spacings  in  each  dimension,  dx  and  dy,  may  differ.  The  elements  are  grouped  hierarchically  at  three  levels. 
The  array  is  subdivided  into  Lx  by  Ly  subapertures,  each  of  which  has  an  associated  time  delay  for 
steering.  Descending  to  the  next  level,  each  subaperture  contains  Mx  by  My  subarrays,  each  of  which 
likewise  has  a  steering  time  delay.  Finally,  each  subarray  contains  Nx  by  Ny  elements,  each  of  which  is 
phase-steered.  Regular  element  spacing  is  maintained  across  subarray  and  subaperture  boundaries.  To 
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array  contains 
L x  by  Ly  time- 
steered  subapertures 


subapertures  contain 
A/x  by  My  time- 
steered  subarrays 


subarrays  contain 
Nx  by  Ny  phase- 
steered  elements 


Figure  3  —  Element  layout  when  diamond  flag  is  nonzero 


simulate  an  array  with  one  time-steered  level  and  one  phase-steered  level,  Lx  and  Ly  should  be  set  to  1;  for 
a  strictly  phased  array,  also  set  Mx  and  My  to  1.  On  the  other  hand,  one  may  model  a  strictly  time-steered 
array  with  one  or  two  hierarchical  levels  by  setting  Nx  and  Ny  to  1  and  dtp  to  360°  (this  renders  the  phase 
shifters  ineffective;  see  Section  4.1). 

Two  additional  program  parameters  specify  special  antenna  geometries.  If  the  flag  diamond  is 
nonzero,  elements  are  effectively  removed  from  even  diagonals  (the  zeroth  diagonal  originates  at  the 
lower-left  element),  leaving  a  diamond  pattern  of  active  elements  as  in  Figure  3.  The  parameter 
azmthOf  f  st  rotates  the  antenna  about  the  z  axis  (boresight);  it  is  the  angle  of  the  positive  x  axis  of  the 
array  above  the  azimuth  reference.  Within  ARRSTATS,  all  calculations  are  performed  in  the  antenna 
coordinate  system;  the  steering  vector  provided  by  the  user  is  transformed  to  the  antenna  coordinate 
system  before  processing,  and  output  coordinates  and  plotted  structures  are  transformed  from  the  antenna 
coordinate  system  after  processing. 
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Finally,  two  parameters  supply  frequency  information.  The  reference  or  design  frequency/,  is  that 
at  which  the  time  and  phase  delays  would  correctly  steer  the  antenna  in  the  absence  of  error  and 
quantization.  The  operating  frequency  /  is  that  at  which  the  simulation  is  to  be  performed.  Multiple 
frequencies  may  be  considered  sequentially  as  described  in  Section  7. 

4.  EXCITATIONS 


4.1.  Error-Free  Excitations 


We  now  construct  the  element  excitations  in  detail  to  show  the  structure  of  the  array  model, 
beginning  with  the  error-free  case  {338-410}.  Because  only  one  frequency  is  considered  at  a  time,  time 
delays  may  be  expressed  as  equivalent  phase  delays.  We  therefore  decompose  the  excitations  into 
magnitude  and  phase  as 


ann,  =\a,,„  lexp(-/0  ),  nwe{0,l,...,LwMwNw-l},  we{x,y}. 


(8) 


where  the  overbars  indicate  the  error-free  case  and  the  nw  label  the  elements  across  the  face  of  the  array, 
ignoring  subaperture  and  subarray  boundaries.  The  error-free  magnitudes  are  made  equal  for  all  active 
elements  and  normalized  to  unit  total  power,  so  that 


(9) 


The  phases  are  derived  from  the  condition  that  at  the  reference  frequency  the  far-field  radiation 
must  interfere  perfectly  constructively  in  the  direction  of  the  steering  vector.  This  implies  that  the  phase 
must  progress  linearly  across  the  face  of  the  array  as 

=  ~ko(dA\  +  dysyny)  +  const.,  (10) 

where  k0  =  2 xf0  !c  is  the  reference  wave  number  and  (sx,  sy)  are  the  direction  cosines  of  the  steering 
vector  {96-166}.  Considering  the  architecture  of  the  array,  the  phases  must  be  built  up  from  the 
equivalent  time  delays  at  the  subaperture  and  subarray  levels  and  the  phase  delays  at  the  element  levels. 
Based  on  Eq.  (10),  the  phase  step  in  direction  w  between  adjacent  elements  in  a  subarray  must  be 

&<Pw  =  ~kod^w-  OD 


Likewise,  since  each  subarray  contains  Nx  by  Ny  elements,  adjacent  subarrays  within  a  subaperture  must 
have  a  phase  difference  of  NyJ^(pw,  which  is  equivalent  to  a  time  step 


Mw=—Nwk<pw 

®o  (12) 

=  --Nwdwsw, 
c 


where  co0  =  27if0  is  the  reference  angular  frequency,  and  the  time  step  across  subapertures  follows 
similarly  as 
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A  Tw  =—NJ4jijw.  (13) 

c 

In  practical  arrays,  the  time  and  phase  delays  are  often  quantized,  leading  to  violations  of  Eq.  (10) 
for  general  steering  angles.  We  suppose  that  the  beamformer  is  capable  of  mitigating  the  effects  of 
subaperture  and  subarray  quantization  by  adjusting  the  subarray  and  element  delays.  That  is,  the  delay 
lost  or  gained  in  each  subaperture  due  to  quantization  can  be  balanced  by  additional  or  reduced  delay  in 
the  subarrays,  assuming  that  the  quantization  interval  of  the  subarrays  is  less  than  that  of  the  subapertures. 
Likewise,  the  error  due  to  subarray  quantization  can  be  balanced  by  adjusting  the  element  phases,  subject 
to  a  similar  assumption.  To  exhibit  this  scheme  mathematically,  we  introduce  the  subaperture,  subarray, 
and  element  quantization  intervals  ST,  St,  and  S<p,  respectively.  We  also  introduce  subaperture  labels  lx  and 
/  and  subarray  labels  mx  and  my;  as  the  nw  ignore  subaperture  and  subarray  boundaries,  so  the  mw  ignore 
subarray  boundaries.  More  specifically,  we  obtain  the  lw  and  mw  from  the  nw  according  to 


/*€ {0,1,.. .1^-1},  and 
mwe{0,l,...LwMw-l}, 

where  |_jcJ  is  the  greatest  integer  not  exceeding  x. 

We  define  the  subaperture  time  delays  without  quantization  or  error  to  be 

=  lxATx  +  lyATy  > 


K  = 

nv 

mxir  = 

tlw 

w 

_NW_ 

(14) 


(15) 


where  the  lw  are  implicitly  dependent  upon  the  nw,  and  the  quantized  subaperture  time  delays  are  then 

T„  n  =  ST  roundif,,  ,,  / ST),  (16) 


where  round  (x)  is  the  integer  nearest  x.  The  subarray  time  delays  contain  an  additional  term  that 
compensates  for  the  subaperture  quantization: 


=  OTX  A*x  +  myAty  —  Tnxlt,  , 

V„  =  #  round  (7/<50, 


(17) 


where  the  mw  are  implicitly  dependent  upon  the  nw,  and  the  element  phase  delays  contain  two  similar 
terms: 


9n,Ky  =  «xA<Px  +  nyA<Py  -  ®0  (tA  +  )>  (]8) 

=S<P  round  ((pli  njs<p) 

In  so  defining  the  delays,  we  have  implicitly  chosen  the  constant  in  Eq.  (10)  to  be  zero.  This  choice 
implies  that  the  lowest  and  leftmost  components  (those  with  lw  -0,mw-  0,  or  nw  =  0)  have  no  associated 
delay  regardless  of  the  steering  vector,  whereas  the  highest  and  rightmost  components  (having  lw  = 
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Lw-  1,  mw  mod  Mw -  M„-  1,  or  mod  AT,,  =  iVw -  1)  have  delays  that  depend  strongly  on  the  steering 
vector. 

In  the  program,  quantization  may  be  avoided  by  setting  the  quantization  steps  to  zero.  The  above 
equations  are  then  equivalent  to 


T„s„y  =i,ly  =ixATx+iyATy, 

tn  n  =  /„  „v  =  (mx  mod  Mx  )Atx  +  (my  mod  My  )Aty ,  and  (19) 

9„x„y  =K»y  =(nxmodNx)A<px+(nymodNy)A<py. 

We  note  that  mw  mod  Mw  is  the  index  of  subarray  mw  within  its  parent  subaperture,  and  likewise 
nw  mod  Nw  is  the  index  of  element  nw  within  its  parent  subarray. 

For  each  element,  the  net  (possibly  quantized)  phase  at  the  operating  frequency  is  the  sum  of 
equivalent  phase  contributions  from  the  three  hierarchical  levels: 


4> 

»X», 


—  vwv  4ywv  ) 


(20) 


where  m  =  2irf  is  the  operating  angular  frequency.  This  quantized  phase  assumes  the  place  of  O  in  Eq. 
(8).  At  the  reference  frequency  (co  =  a>0 )  and  with  no  phase  quantization,  this  construction  of  the  phase 
yields  the  linear  progression  of  Eq.  (10). 


4.2.  Erroneous  Excitations 

We  model  the  errors  in  a  real  antenna  by  applying  amplitude  and  phase  errors  to  each  level  of  the 
antenna  hierarchy  {525-549}.  Amplitude  errors  multiply  the  error-free  amplitudes  by  factors  of  the  form 
1  +p  where  p  is  a  random  number,  distributed  normally  with  zero  mean.  Each  level  contributes  such 
errors,  so  that  the  erroneous  amplitude  for  element  (nx,  ny)  is 


I  anjK I  _  I  a»,n  I  0  +  I  )0  +  +  Phxh.  )  » 


(21) 


where  R,  r,  and  p  correspond  to  subapertures,  subarrays,  and  elements,  respectively.  This  model  allows 
corresponding  elements  in  different  subarrays  to  contribute  distinct  errors,  and  likewise  for  corresponding 
subarrays  within  different  subapertures.  The  user  specifies  the  standard  deviations  <rR,  <rr,  and  crp  of  the 
respective  amplitude  errors. 

Time  and  phase  errors  add  to  the  error-free  (but  possibly  quantized)  time  and  phase  delays.  The 
subapertures  and  subarrays  contribute  random  time  errors  T  and  t,  respectively,  and  the  elements 
contribute  random  phase  errors  (p,  all  drawn  from  zero-mean  normal  distributions.  The  user  specifies  the 
corresponding  standard  deviations  <rT,  <r„  and  o^.  Additionally,  the  user  may  specify  a  deterministic  time 
error  .T  for  each  subaperture.  The  net  equivalent  phase  error  for  element  («x,  ny)  is 


(22) 


and  the  total  erroneous  phase  is  the  sum  of  the  quantized  and  error  phases: 
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(T)  =<t>  -f 

"A  «xWy 


(23) 


Finally,  the  erroneous  complex  excitations  are 

exP(-/(I>«x»y  >  ■ 


(24) 


5.  RADIATION  PATTERN 

In  standard  array  theory  with  mutual  coupling  ignored,  the  field  pattern  is  the  product  of  the  array 
factor  and  the  element  factor.  In  ARRSTATS,  the  element  factor  is  unity,  corresponding  to  uniform 
hemispherical  radiation,  so  the  field  pattern  equals  the  array  factor.  (See  the  code  {2595-2715}  for  notes 
on  expanding  ARRSTATS.)  Given  the  complex  excitations,  the  array  factor  (and  field  pattern)  in  the 
direction  (£  17)  is  given  by  the  two-dimensional  Fourier  transform  {551-561} 


g(&i/)  =  X IexP[-'^"x  +tjdyny )]a/lx„y  , 

»x  "y 


(25) 


where  k  =  lirflc  is  the  operating  wave  number.  ARRSTATS  uses  a  fast  Fourier  transform  to  obtain  g  in. 
discrete  directions  given  by 


K0x’ 


MyQy’ 


?xe{0,l,...,ex-l},  and 

qye{0,l,...,Qy -1}, 


(26) 


where  Qx  and  Q  are  the  number  of  points  in  the  transform  in  each  dimension  {168-173}.  The  (£,x»fyy) 
grid  is  extrapolated  to  all  of  visible  space  using  {412-497,  558} 


+&,’  7(/y  )  £(‘=<7,  ’  )  and 

+Qy  )  _  5  )> 


which  are  valid  for  all  integers  qx  and  qy.  Because  of  the  normalization  condition  of  Eq.  (9), 

|g(£»/)l^1  ^ 

for  all  £  and  tj,  with  the  equality  holding  only  where  the  excitations  interfere  perfectly  constructively. 
Therefore,  when  the  field  pattern  is  expressed  in  decibels,  0  dB  corresponds  to  perfectly  constructive 
interference. 

6.  PATTERN  MEASURES 

The  code  at  the  heart  of  the  program  analyzes  the  radiation  pattern  to  obtain  several  measures  that 
quantify  the  characteristics  and  performance  of  the  array.  The  following  subsections  describe  these 
measures,  generally  focusing  on  the  concepts  behind  them  and  on  their  interpretation  rather  than  on  the 
specific  method  of  calculation.  Details  of  the  algorithms  and  further  discussion  may  be  found  in  the  code. 
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6.1.  Pointing  Vector,  Pointing  Error,  and  Peak  Power  Density 

One  of  the  most  significant  pattern  measures  is  the  direction  of  maximum  radiation,  here  called 
the  pointing  vector.  The  program  determines  it  directly  from  the  field  pattern  and  also  indirectly  from  the 
excitations;  the  final  pointing  vector  is  a  weighted  average  of  the  two,  as  discussed  below.  In  identifying 
the  pointing  vector  from  the  field  pattern  {563-787},  the  program  first  locates  the  pattern’s  maximum 
magnitude  over  the  grid  of  discrete  direction  cosines.  For  a  well-formed  beam,  the  neighboring  samples 
should  fall  off  parabolical ly,  so  they  are  fitted  to  the  elliptic  paraboloid13 

|g(£*)l  +W^  +  \vtl1  +X£  +  Ytj  +  Z  (29) 

in  a  least-squares  sense.  If  XJV >  IT3,  as  should  be  the  case  for  a  normal  beam,  the  conic  is  indeed  elliptical 
(corresponding  to  its  level  curves),  and  its  maximum  occurs  at  (<f,  tj)  =  (plx,  ply),  where  pix  and  ply  satisfy 


U  wY. Plx 
w  vlPxy 


(30) 


The  coordinate  pair  (plx,ply)  is  taken  to  be  the  pointing  vector  for  the  first  method.  The  deviation  of  the 
pattern  samples  from  conic  form  is  used  to  construct  a  covariance  ellipse  for/7Ix  and  ply  that  expresses  the 
uncertainty  in  their  values. 

The  second  method  determines  the  pointing  vector  from  the  complex  excitations  {789-961}. 
Because  a  well-behaved  array  will  have  a  nearly  flat  phase  front,  the  excitation  phases  On  ^  are  fit  to  the 
plane  -knxdj  -  knydyr\  -  d>0,  weighted  according  to  the  excitation  magnitudes.  The  pointing  vector 
coordinates  (p2x,/»2y)  are  the  direction  cosines  (<f,  rj)  that  give  the  best  fit.  As  with  the  first  method,  a 
covariance  ellipse  for  p2x  and  p2y  is  obtained  from  a  measure  of  the  deviation  from  the  plane. 

Each  method  is  useful  but  limited.  The  first,  the  fit  of  the  transform,  is  robust  even  for  poorly 
aimed  and  malformed  beams,  but  it  is  limited  by  the  resolution  of  the  Fourier  transform.  The  second,  the 
fit  of  the  excitations,  is  independent  of  transform  resolution  but  accurate  only  for  nearly  planar  phase 
fronts,  approaching  the  exact  solution  as  the  phase  and  time  errors  and  quantization  intervals  decrease.  To 
obtain  a  single  pointing  vector  (px,py),  the  program  averages  pointing  vectors  from  the  two  methods, 
weighting  each  by  the  inverse  of  the  area  of  its  covariance  ellipse  {963-990}.  An  average  covariance 
ellipse  is  also  constructed  in  a  consistent  manner. 

With  the  final  pointing  vector  in  hand,  the  pointing  error  y  is  straightforwardly  obtained  {992- 
1046}  from 


cosy  =  s-p 

=  sxpx  +  sypy  +  szpz. 


(31) 


where  sz  and  pz  follow  from  unit  magnitude  constraints  on  s  and  p.  If  one  desires  the  uncertainty  in  y  due 
to  uncertainty  in  p,  an  alternate  calculation  based  on 

siny  =  |sx  p|  (32) 


may  be  selected  in  the  program. 
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The  peak  power  density 


Pma*=\g(Ps’Py)\2 


(33) 


is  calculated  directly  from  Eq.  (25). 

6.2.  Main  Beam  Region 

The  main  beam  region  is  the  set  of  field  pattern  samples  that  are  judged  to  belong  to  the  main 
beam.  Although  normally  not  of  interest  as  a  final  measure,  it  is  essential  for  obtaining  subsequent 
measures.  It  may  be  constructed  conceptually  by  imagining  a  contour  at  an  adjustable  level.  Beginning  at 
the  pattern  maximum,  we  decrease  the  level  so  that  the  contour  expands  in  size,  following  the  topology  of 
the  main  beam.  Eventually  the  contour  will  intersect  a  local  minimum  or  a  saddle  point;  the  closed 
contour  about  the  pattern  maximum  at  that  level  delineates  the  main  beam  region  inside  from  the  sidelobe 
region  outside.  Equivalently,  that  contour  is  the  lowest  one  containing  the  global  maximum  that  encircles 
no  other  local  maxima.  The  power  level  of  the  contour  is  called  the  beam  depth.  Generally,  well-formed 
beams  are  deep  (that  is,  the  beam  depth  is  much  less  than  one),  while  malformed  beams  are  shallow,  but 
the  user  should  keep  in  mind  that  the  beam  depth  depends  on  the  transform  resolution.  In  the  program 
{1048-1137},  the  beam  depth  is  stored  in  the  variable  beamDepthDB  in  decibels  relative  to  Pmm  and  is 
output  to  the  user.  Information  obtained  while  determining  the  main  beam  region  is  used  in  finding  the 
main  beam  width  and  roll,  below,  and  the  region  itself  is  used  directly  in  calculating  the  powers  radiated 
into  the  main  beam  and  sidelobes. 

6.3.  Main  Beam  Width  and  Roll 

The  level  contours  of  the  main  beam  generated  by  a  two-dimensional  array  are  nominally 
ellipses;14’15  therefore  they  can  be  described  by  their  center  location,  major  and  minor  axes,  and 
orientation.  Having  already  obtained  a  measure  of  the  beam’s  location  in  the  form  of  the  pointing  vector, 
we  use  the  major  and  minor  axes  and  orientation  of  the  elliptical  contour  at  a  given  power  level  to 
describe  the  beam’s  shape  {1139-1403}.  The  conventional  power  level  is  Pmax  / 2,  or  about  3  dB  down. 
The  angle  subtended  by  the  ellipse’s  longest  diameter  —  its  major  axis is  taken  as  the  beam’s  major 
width  (that  is,  full  width  at  half  maximum  power);  that  subtended  by  its  shortest  diameter,  the  beam’s 
minor  width.  The  ellipse  orientation  gives  the  roll  angle,  but  we  must  choose  an  origin  for  the  orientation 
angle.  Construct  three  great  circles  as  in  Figure  4:  A,  along  the  azimuth  reference;  B,  connecting  boresight 
and  the  pointing  vector;  and  C,  along  the  beam’s  major  width.  The  roll  angle  is  defined  as  the  sum  of  the 
angle  between  A  and  B  and  that  between  B  and  C.  It  happens  that  the  roll  angle  so  defined  is  merely  the 
orientation  of  the  major  width  relative  to  the  azimuth  reference  when  viewed  in  the  stereographic 
projection,  which  preserves  angles  between  great  circles.  Moreover,  because  the  great  circles  along  the 
beam’s  major  and  minor  widths  intersect  orthogonally  on  the  hemisphere,  their  stereographic  projections 
do  also.  These  facts  motivate  the  program’s  use  of  stereographic  coordinates  for  fitting  an  ellipse  to  the 
level  contour  and  determining  its  major  and  minor  axes  and  orientation.  However,  the  roll  angle  and  beam 
widths  thereby  obtained  are  approximate  for  two  reasons.  First,  for  beams  off  boresight,  the  great  circles 
along  the  major  and  minor  widths  project  as  circles,  whereas  the  major  and  minor  axes  of  the  projected 
ellipse  are  straight  line  segments.  Second,  the  local  scale  in  the  projection  increases  away  from  boresight, 
artificially  enlarging  the  portion  of  the  beam  farthest  from  boresight.12  The  error  may  become  significant 
for  beams  far  from  boresight.  A  more  sophisticated  method  of  determining  the  beam  widths  is  suggested 
near  the  end  of  the  program  code. 


12 


C.  S.  West 


Figure  4  —  Beam  width  ellipse  and  roll  angle.  The  great  circles  A,  B,  and  C  are  described  in  the  text,  and  the  angles  between 
them,  indicated  with  thick  arcs,  are  added  to  obtain  the  roll  angle.  The  left  projection  is  orthographic;  the  right,  stereographic,  i 
The  spherical  coordinates  of  the  pointing  vector  are  0  =  40°  and  y/  =  20°,  the  beam’s  major  and  minor  widths  are  16°  and  8°,  and  , 
the  roll  angle  is  90°.  r 


6.4.  Directivity,  Power  Ratios,  and  Average  Sidelobe  Level 

The  next  measures  obtained  all  depend  on  integrals  of  the  power  pattern  over  solid  angle  regions 
{1405-1527}.  The  integrals  are  calculated  from  the  discrete  samples  of  the  pattern  using  the  midpoint 
approximation,  as  detailed  in  the  code.  The  solid  angle  regions  of  interest  are  visible  space  (the  half-space 
z  >  0);  the  main  beam,  as  described  by  the  main  beam  region,  above;  and  the  sidelobes,  defined  as  all 
regions  of  visible  space  not  in  the  main  beam.  The  program  determines  the  total  powers  radiated  into 
these  regions;  call  them  for  visible  space,  n„,  for  the  main  beam,  and  IIS  for  the  side  lobes.  We  have 

nv=nm  +ns.  (34) 

The  directivity  is  the  ratio  of  maximum  to  average  power  densities,  assuming  no  back  radiation  into  z<  0: 

[)  -  ^max  (35) 

IIv/4jr’ 


The  program  also  calculates  the  power  ratios  n,n  /  IIV,  nm  /ns,  and  riv  /Fls,  which  generally  decrease  as 
the  beam  degrades.  Finally,  the  average  sidelobe  level  relative  to  the  beam  peak  is 


L 

avE  a 


(36) 


where  Qs  is  the  solid  angle  occupied  by  the  sidelobes. 
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6.5.  Powers  and  Distances  of  Largest  and  Nearest  Sidelobes 

The  last  analysis  identifies  the  sidelobe  with  the  largest  power  density  and  the  sidelobe  closest  to 
the  main  beam,  and  for  each  it  reports  the  power  density  and  angular  distance  from  the  main  beam 
{1529-1677}.  The  power  levels  and  locations  of  the  sidelobes  are  obtained  by  fitting  to  elliptic 
paraboloids  (Eq.  (29)),  and  the  power  levels  are  reported  relative  to  the  beam  peak  Pmax. 

7.  MULTIPLE  REALIZATIONS  AND  PARAMETER  VALUES 

The  analyses  described  above  apply  to  the  radiation  pattern  corresponding  to  a  single  set  of  array 
parameters  and  one  realization  of  any  random  variables.  ARRSTATS  contains  two  outer  loops  with  which  it 
analyzes  multiple  radiation  patterns;  one  is  a  loop  over  realizations  of  random  variables,  and  the  other 
loops  over  a  user-selected  independent  variable.  The  loop  over  realizations  {190-218,  521-523,  1679- 
1819, 1905}  is  motivated  by  the  following:  When  simulating  random  errors,  one  is  usually  interested  not 
in  the  performance  obtained  by  one  realization  of  the  errors  but  rather  in  the  performance  statistics  for  an 
ensemble  of  statistically  identical  arrays.  To  that  end,  ARRSTATS  can  generate  an  arbitrary,  user-specified 
number  of  realizations  for  which  it  will  accumulate  performance  statistics.  The  program  outputs  the  mean 
and  standard  deviation  of  each  performance  measure.  The  second  loop  {175-188,  326-336,  2578-2579}, 
over  an  independent  variable,  allows  the  user  to  examine  the  variation  of  performance  measures  as  the 
variable  changes.  Possible  independent  variables  include  but  are  not  limited  to  the  operating  and  reference 
frequencies,  steering  angles,  error  parameters,  quantization  intervals,  and  even  parameters  of  the  array 
geometry.  ARRSTATS  produces  a  graph  showing  each  measure  as  a  function  of  the  independent  variable, 
as  discussed  below. 

8.  PROGRAM  OUTPUT 

ARRSTATS  outputs  its  results  in  three  ways:  textual  output  of  running  statistics,  a  plot  of  the 
radiation  pattern  for  the  last  realization  in  an  ensemble,  and  a  summary  plot  of  the  performance  measures 
as  functions  of  the  independent  variable.  The  textual  output  is  a  table  like  that  in  Fig.  5  printed  to  the 
MATLAB  command  window  {1864-1903}.  The  values  in  the  table  are  the  means  and  standard  deviations 
of  the  performance  measures  for  all  members  of  the  statistical  ensemble  that  have  been  realized  thus  far. 
The  table  may  be  interpreted  according  to  the  descriptions  in  Section  6,  keeping  in  mind  the  following. 


Means  and  [std  devs]  for  16  of  16  realizations 


beam  direction 
pointing  error 
peak  power  dens 
beam  depth 
beam  width 
beam  roll 
directivity 
power  ratio  m/v 
power  ratio  m/s 
power  ratio  v/s 
avg  sidelobe 
nrst  sidelobe 
lgst  sidelobe 


(29.954,  60.007)  deg,  std  dev  0.038  deg 


0.0528 

[0.0279)  deg 

-0.426 

[0.034 

]  dB 

-24.66 

[1.49 

]  dB  re  peak 

(  2.095 

[0.005 

],  1.813  [0.004]) 

deg 

59.56 

[0.65 

]  deg 

38.894 

[0.034 

]  dB 

-1.360 

[0.030 

]  dB 

4.347 

[0.113 

)  dB 

5.707 

[0.083 

3  dB 

-41.61 

[0.11 

]  dB  re  peak 

-13.16 

[0.69] 

dB  re  peak,  3.00 

[0.01] 

deg 

off 

beam 

-12.30 

[0.58] 

dB  re  peak,  3.11 

[0.09] 

deg 

off 

beam 

Figure  5 _ Textual  output  of  running  statistics.  The  parameters  of  this  example  are  in  the  code  listing;  for  the  values  above,  the 

standard  deviation  of  the  subarray  time  error  is  10  ps  (stdTimeMPS  =  10). 
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The  coordinates  specifying  the  beam  direction  are  the  spherical  coordinates  ( 6  and  y/)  of  the  average 
pointing  vector  {1821-1862}.  The  standard  deviation  of  the  beam  direction  is  an  estimate  of  the  rms 
angular  deviation  of  the  ensemble  of  pointing  vectors  from  their  mean  {1821-1862}.  Both  the  beam 
direction  coordinates  and  roll  angle  include  compensation  for  the  azimuth  offset  azmthOf  f  st.  The  peak 
power  density  is  relative  to  perfectly  constructive  interference,  and  the  beam  depth,  average  sidelobe 
level,  and  levels  of  the  nearest  and  largest  sidelobes  are  relative  to  the  peak  power  density  (suggested  by 
the  use  of  “dB  re  peak”  in  the  table). 

To  aid  in  visualizing  an  array’s  performance,  ARRSTATS  can  graphically  present  the  radiation 
pattern  and  several  of  the  performance  measures  as  in  Fig.  6  {220-324,  2138-2574}.  A  significant 


Figure  6  —  Lambert  projection  of  the  radiation  pattern  with  a  spherical  coordinate  grid  superimposed.  Boresight  is  at  the  center, 
and  the  steering  vector  is  (0  =  30°,  y/  =  60°).  The  white  dot  denotes  the  beam  peak;  *,  the  largest  sidelobe;  and  +,  the  nearest 
sidelobe.  The  legend  gives  the  power  in  dB  relative  to  perfectly  constructive  interference.  The  pattern  is  one  realization  of  the 
case  stdTimeMPS  -  10  for  the  parameters  given  in  the  code  listing. 
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plotting  option  is  the  choice  of  projection  from  among  those  described  in  Section  2.  Briefly,  the  Lambert 
projection  preserves  the  relative  areas  of  regions  (for  example,  the  size  of  the  main  beam  relative  to  the 
hemisphere  or  to  prominent  sidelobes),  the  stereographic  projection  preserves  local  shape  (and  the 
orthogonality  of  the  major  and  minor  beam  width  cuts),  and  the  orthographic  projection  gives  a  picture  of 
the  hemisphere.  The  pattern  may  be  plotted  linearly  or  logarithmically  in  power,  and  the  logarithmic 
depth  may  be  specified.  The  user  may  also  specify  the  color  map  and  shading  method  to  use.  The 
performance  measures  that  can  be  graphically  indicated  on  the  radiation  pattern  include  the  pointing 
vector,  the  axes  of  the  uncertainty  ellipse  of  the  pointing  vector,  the  actual  and  fitted  beam  width  contours, 
the  main  beam  region,  the  locations  of  the  nearest  and  largest  sidelobes,  and  other  measures  of  less 
frequent  interest.  If  multiple  realizations  are  generated,  the  radiation  pattern  will  be  plotted  only  for  the 
last  realization  as  a  representative  of  the  ensemble. 

If  the  outer  loop  over  an  independent  variable  is  used,  the  summary  figure  plots  the  performance 
measures  as  functions  of  the  independent  variable,  as  in  Fig.  7  (220-324,  1907-2136}.  The  figure  groups 
fifteen  measures  (all  except  the  beam  direction)  into  eight  subplots  and  utilizes  distinct  colors  or  line 
styles  and  both  left  and  right  axes  for  the  ordinates.  The  title  of  each  subplot  gives  the  names  of  the 
measures;  where  multiple  colors  or  linestyles  appear,  each  measure  name  is  followed  by  a  color  or  style 
code  in  parenthesis,  and  where  left  and  right  axes  are  used,  an  axis  code  (“L”  or  “R”)  also  appears.  If 
more  than  one  realization  was  generated  for  each  value  of  the  independent  variable,  error  bars  extend  one 
standard  deviation  above  and  one  standard  deviation  below  each  point. 

9.  SUMMARY 

We  have  introduced  a  computer  program  for  assessing  the  effects  of  errors  in  rectangular-grid 
planar  array  antennas.  The  most  general  array  comprises  phase-steered  elements  grouped  into  time- 
steered  subapertures  and  subarrays;  this  includes  strictly  phase-steered  arrays  and  time-steered  arrays  as 
particular  cases.  The  time  and  phase  delays  may  be  quantized,  and  random  errors  may  be  assigned  to  the 
times,  phases,  and  amplitudes.  From  simulated  radiation  patterns,  the  program  obtains  performance 
measures  over  statistical  ensembles  and  as  functions  of  a  user-specified  independent  variable,  producing 
textual  and  graphical  output. 


Appendix  A 

LIST  OF  VARIABLES 


This  list  of  significant  variables  mentioned  in  this  report  gives  their  symbol  in  this  report,  coded 
name,  and  description.  It  does  not  include  all  program  variables.  Ellipses  (...)  stand  for  prefixes,  and 
asterisks  (*)  indicate  that  the  variable  holds  the  specified  quantity  temporarily. 


Text 

Code 

Description 

x,y,z 

coordinates  in  real  space 

dirCosX,  ...Y,  ...Z 

direction  cosines 

SX,  Sy 

sx,  sy 

steering  vector  direction  cosines 

e 

...Polar 

spherical  polar  coordinate 

¥ 

...Azmth 

spherical  azimuth  coordinate 

a 

..Az 

traditional  azimuth  coordinate 

e 

...El 

traditional  elevation  coordinate 

r 

radial  coordinate  of  projection 

dy 

dx,  dy 

element  spacings 

numLX,  numLY 

numbers  of  subapertures 

My 

numMX,  numMY 

numbers  of  subarrays  per  subaperture 

^Ny 

numNX,  numNY 

numbers  of  elements  per  subarray 

4’  4 

lx,  ly 

subaperture  labels 

mx,  my 

mx,  my 

subarray  labels 

nx,riy 

nx,  ny 

element  labels 

diamond 

indicates  diamond  element  pattern 

azmthOf fst 

antenna  rotation  angle  about  boresight 

fo 

fRef 

reference  frequency 

f 

fOpr 

operating  frequency 

K 

reference  wave  number 

k 

operating  wave  number 

co0 

reference  angular  frequency 

CO 

operating  angular  frequency 

stdAmplL 

standard  deviation  of  subaperture  amplitude  error 

stdAmplM 

standard  deviation  of  subarray  amplitude  error 

Op 

stdAmplN 

standard  deviation  of  element  amplitude  error 

R 

random  subaperture  amplitude  error 

r 

random  subarray  amplitude  error 

17 


18 


C.  S.  West 


Text 

Code 

P 

a 

excMagldl 

a 

excMag 

ST 

qntTimeL 

St 

qntTimeM 

Stp 

at;,  at; 

A/x,  A ty 

A <PX,  A <py 

qntPhseN 

T 

"x"y 

TimeL* 

TimeM* 

^Mx/Jy 

PhseN* 

T 

»xny 

TimeL 

tfixny 

TimeM 

*Pnxny 

PhseN 

stdTimeL 

stdTimeM 

*<p 

stdPhseN 

jr 

of sTimeL 

f 

7 

9 

O 

0> 

excPhsIdl 

$ 

excPhsErr 

<D 

excPhs 

g 

g 

Q»Qy 

tx,  ty 

Px’Py 

px,  py 

y 

errPoint 

p 

max 

gSqrMax 

Aivg 

powerSideAvgDB 

nv 

powerVisb 

nin 

powerMain 

ns 

powerSide 

D 

directivityDB 

Description _ 

random  element  amplitude  error 
error-free  (ideal)  excitation  magnitude 
erroneous  (actual)  excitation  magnitude 
subaperture  time  quantization  interval 
subarray  time  quantization  interval 
element  phase  quantization  interval 
subaperture  time  steps 
subarray  time  steps 
element  phase  steps 
error-free  subaperture  time  delays 
error-free  subarray  time  delays 
error-free  element  phase  delays 
quantized  subaperture  time  delays 
quantized  subarray  time  delays 
quantized  element  phase  delays 
standard  deviation  of  subaperture  time  error 
standard  deviation  of  subarray  time  error 
standard  deviation  of  element  phase  error 
deterministic  subaperture  time  error 
random  subaperture  time  error 
random  subarray  time  error 
random  element  phase  error 
error-free  excitation  phase 
quantized  excitation  phase 
excitation  phase  error 
erroneous  (actual)  excitation  phase 
field  pattern 

number  of  points  in  Fourier  transform 
pointing  vector  direction  cosines 
pointing  error 
maximum  power  density 

average  relative  sidelobe  level  (powerSideAvgDB  =  10  log,0  Zavg) 

power  radiated  into  visible  space 

power  radiated  into  the  main  beam 

power  radiated  into  the  side  lobes 

directivity  (directivityDB  =  10  log10  D) 


Appendix  B 

PROGRAM  LISTING 


1  A  Calculate  performance  parameters  of  an  array  with  excitation  errors 

2  A 

3  A  The  excitation  time  and  phase  convention  is  exp  (-i  (omega  t  +  phi)); 

4  A  positive  (negative)  phases  phi  correspond  to  a  leading  (lagging) 

5  V.  excitation.  Distances  are  stored  in  meters;  times,  nanoseconds; 

6  A  frequencies,  gigahertz.  Angles  (both  geometric  and  phase)  are  always 

7  V.  specified  in  radians.  Generally,  the  coordinate  x  increases  with  the 

8  A  column  index;  y,  with  the  row  index. 

9  A 

10  'A  Parameters  that  may  be  changed  by  the  user  are  marked  "INPUT"  in 

11  'A  comments.  Frequently-used  inputs  appear  near  the  top  of  the  program, 

12  A  but  some  inputs  are  defined  elsewhere,  particularly  in  the  plotting 

13  'A  section. 

1 4 

15  V,  Suggestions  for  improvements  are  provided  at  the  end  of  the  program. 

16  A 

17  V.  This  program  is  coded  for  Matlab  version  5.3  (Release  11),  although 

18  A  most  of  the  code  will  run  under  version  4.2. 

19  A 

20  A  Written  by  Stan  West,  1998,  1999 

21  V.  U.S.  Government  work  not  subject  to  copyright 

22 

23  ‘A  Declare  physical  and  conversion  constants 

24  A 

25  c  “  0.299792458;  'A  speed  of  light  in  m/ns 

26  rpd  »  pi  /  180;  A  radians  per  degree 

27  twopi  -  2  *  pi; 

28 

29  'A  Set  operating  frequency  relative  to  reference  frequency 

30 

31  'A  At  the  reference  frequency,  or  fOpr  -  fRef,  the  time  and  phase  delays 

32  ‘A  will  properly  steer  the  antenna  in  the  absence  of  error  and 

33  ‘A  quantization. 

34 

35  fRef  -  3.0;  A  INPUT  reference  frequency  in  GHz 

36  fOpr  -  fRef  +1*0. 5/2;  A  INPUT  operating  frequency  in  GHz 

37  'A  center  frequency  +  (-1  ...  1}  *  bandwidth  /  2 

38 

39  A  Describe  array  geometry  and  error  statistics 

40 

41  A  Elements  lie  on  a  regular  planar  grid  and  are  grouped  heirarchically . 

42  A  The  array  is  subdivided  into  numLX  subapertures  in  the  x  dimension  and 

43  ‘A  numLY  subapertures  in  the  y  dimension.  Each  subaperture  is 

44  A  time-steered  and  has  associated  with  it  a  user-set  deterministic 

45  A  absolute  time  error,  a  random  absolute  time  error,  and  a  random 

46  A  relative  amplitude  error.  The  random  errors  are  normally-distributed 

47  A  with  zero  mean  and  user-set  standard  deviation.  Each  subaperture 

48  A  contains  numMX  by  numMY  subarrays,  each  of  which,  like  the 

49  A  subapertures,  is  time-steered  and  has  a  random  absolute  time  error  and 

50  A  a  random  relative  amplitude  error.  Finally,  each  subarray  contains 

51  A  numNX  by  numNY  elements,  each  of  which  is  phase-steered  and  has  a 

52  A  random  absolute  phase  error  and  a  random  relative  amplitude  error. 

53  A  All  elements  are  spaced  by  dx  in  x  and  dy  in  y,  even  across  subarrays 

54  'A  and  subapertures. 

55 

56  A  Set  parameters  of  subapertures 

57 

58  numLX  =  1;  A  INPUT  number  of  subapertures  in  x 

59  numLY  =  1;  A  and  y  dimensions 

60  qntTimeL  =  0;  A  INPUT  quantization  interval  in  ns;  0  for  no  quantization 

61  ofsTimeL  =  zeros  (numLY,  numLX);  A  INPUT  deterministic  absolute  time  error  in  each  subaperture  in  ns 

62  stdTimeL  ■=  0.000;  |  INPUT  standard  deviation  of  absolute  time  error  in  each  subaperture  in  ns 

63  stdAmplL  =  0.0;  A  INPUT  standard  deviation  of  relative  amplitude  error 

64  A 

65  v.  Set  parameters  of  subarrays 

66 

67  numMX  =  8; 

68  numMY  -  8; 

69  qntTimeM  *  0; 

70  stdTimeM  »  0.00; 

71  StdAmplM  -  0.0; 

72 

73  A  Set  parameters  of  elements 

74  A 

75  dx  *  1.6  *  0.0254;  A  INPUT  element  spacing  in  x 


INPUT  number  of  subarrays  in  x 
A  and  y  dimensions  per  subaperture 

A  INPUT  quantization  interval  in  ns;  0  for  no  quantization 
A  INPUT  standard  deviation  of  absolute  time  error  in  each  subarray  in  ns 
A  INPUT  standard  deviation  of  relative  amplitude  error 
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76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 
167 


dy  ■=  1.6  *  0.0254; 
numNX  =  8 ; 

numNY  =  8 ; 

qntPhseN  =  0  *  rpd; 

stdPhseN  =  0  *  rpd; 

stdAmplN  =  0.0; 


and  y  dimensions  in  meters 
INPUT  number  of  elements  in  x 
and  y  dimensions  per  subarray 

INPUT  quantization  interval  in  radians;  0  for  no  quantization 

INPUT  standard  deviation  of  absolute  phase  error  in  each  element  in  radians 

INPUT  standard  deviation  of  relative  amplitude  error 


Set  other  parameters 

diamond  =  0;  INPUT  0:  full  array;  1:  eliminate  excitations  on  even  diagonals 

azmthOf fst  =  0  *  rpd;  v.  INPUT  angle  of  the  positive  x  axis  of  the  array  above  the  azimuth  reference 


V.  Alternatively/  select  an  array 


switch  0  INPUT  case  number  for  arrays  below  or  0  to  use  values  above 
case  1 

!?.  insert  frequency,  array,  and  error  parameters  here 
case  2 
end 


V.  Specify  steering  angle 
V. 

V.  Two  coordinate  systems,  described  below,  are  available  for  specifying 
V,  the  steering  angle:  traditional  azimuth/elevation  coordinates  and 

spherical  coordinates.  Azimuth/elevation  coordinates  are  converted  to 
%  spherical  coordinates  for  internal  program  use,  and  output  is  given  in 
%  spherical  coordinates. 

% 

V.  In  the  diagrams  below,  the  antenna  lies  in  the  x*-y*  plane  with 
V.  boresight  along  the  positive  z*  axis.  The  antenna's  z  axis  coincides 
V,  with  z  * ,  and  its  x  axis  is  at  an  angle  azmthOf  fst  above  the  x'  axis, 

V,  which  is  the  azimuth  reference. 

*  The  equations  relating  the  the  Cartesian  coordinates  x’,  y’,  and  z'  of 
%  a  unit  vector,  the  spherical  coordinates  Polar  and  Azmth,  and  the 
?.  traditional  coordinates  Az  and  El  are 
V. 

•t  x*  *  sin  Polar  cos  Azmth  *  -cos  El  sin  Az 

y»  -  sin  Polar  sin  Azmth  -  sin  El 
V,  z*  -  cos  Polar  -  cos  El  cos  Az 


?.  cos  Polar  *=  cos  El  cos  Az  »  z' 

V.  tan  Azmth  -  -tan  El  /  sin  Az  «  y*  /  x* 


% 

’A  -tan  Az  =  tan  Polar  cos  Azmth  - 

V,  sin  El  “  sin  Polar  sin  Azmth  - 

% 

%  Since  (x',  y',  z*)  is  a  unit  vector, 
%  cosines. 

% 

switch  2 
case  1 

4  y* 

v,  I 

V.  —  I 

¥.  /  I 

main  _  /  I 

V.  beam  — _  I 

%  /  —  o - 

El  I  /  x’ 

}_-  / 

V. 

¥.  proj .  in  / 

'l  xz  plane  Az  z* 

steerAz  **  30  *  rpd; 

steerEl  =  30  *  rpd; 

steerPolar  -  acos  (cos  (steerEl) 
steer Azmth  »  atan2  (tan  (steerEl) 


x*  /  z  * 

y* 

its  components  are  direction 


Use  azimuth/elevation 
coordinates.  Given  the 
projection  of  the  steering  vector 
onto  the  x*-zf  plane,  the  azimuth 
is  the  angle  between  it  and  the 
z'  axis,  while  the  elevation  is 
the  angle  between  it  and  the 
steering  vector.  In  the  diagram, 
both  angles  are  positive. 


?,  INPUT  azimuth  angle 
V.  INPUT  elevation  angle 
r  cos  (steerAz)); 

-sin  (steerAz)); 


case  2 

if,  proj  -in  y ' 

V,  xy  plane  I 

\_—  I 

\  I  \  Azmth 

V.  main  _  /\  !  \ 

V.  beam  --/_  \l  I 

I  I  — o- 
I  I  / 

V.  Polar  \  \/ 

\/ 

/ 


steerPolar  -  30  *  rpd; 
steerAzmth  =  60  *  rpd; 

otherwise 

error  {'Invalid  switch  parameter.') 

end 

steerAzmth  =  steerAzmth  -  azmthOffst; 


Use  spherical  coordinates.  The 
polar  angle  is  the  angle 
between  the  z'  axis  and  the 
steering  vector.  The  azimuth 
angle  is  the  angle  between  the 
x*  axis  and  the  projection  of 
the  steering  vector  onto  the 
x'-y’  plane.  In  the  diagram, 
both  angles  are  positive. 


‘i.  INPUT  polar  angle 
V.  INPUT  azimuth  angle 


V.  now  relative  to  antenna's  x  axis 
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237 
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7.  Set  transform  size 

l  (see  comments  elsewhere  related  to  the  discrete  Fourier  transform.) 


tx 

ty 


2A9; 

2A9; 


INPUT  number  of  transform  points  in  x 
and  y 


7.  Declare  independent  variable  and  its  values 


The  main  loop  iterates  over  values  of  the  independent  variable  named 
below  The  independent  variable  may  be  any  parameter,  including,  for 
example,  those  describing  geometry,  frequency,  error,  quantization, 
and  steering  angle.  It  may  also  be  an  otherwise  unknown  variable  that 
is  transformed  to  a  known  program  parameter  by  custom  code  m  the  main 
,  loop.  To  effectively  disable  the  loop,  set  a  dummy  variable  to  a 
,  scalar  value. 


V.  INPUT  name  of  independent  variable 
7.  INPUT  vector  of  values  it  will  assume 
7  make  it  a  column  vector 


V. 

indVarName  -  'stdTimeMPS* ; 
indVar  -  0  :  2  :  20; 
indVar  -  indVar  (:); 
indVarLen  -  length  (indVar); 

V,  Initialize  statistics  variables 

%  If  the  excitations  are  random,  one  is  often  interested  in  the  mean  and 
V.  standard  deviation  of  the  performance  measures.  For  each  value  of  the 
independent  variable,  the  program  will  generate  numRlz  realizations  of 
Vi  the  random  excitations  and  accumulate  the  statistics  of  the 
V.  performance  measures. 


numRlz  *=  16; 

numAcc 

pxS 

pyS 

pzS 

errPointS 

V.  errPointUncS 

beamPowerDBS 

beamDepthDBS 

hpbwMjrS 

hpbwMnrS 

rolls 

di recti vi tyDBS 

powerMainVisbDBS 

powerMainSideDBS 

powerVisbSideDBS 

powerSideAvgDBS 

slNrstDistS 

slNrstPowrDBS 

slLgstDistS 

slLgstPowrDBS 


nan  *  ones  (indVarLen,  1); 
nan  *  ones  (indVarLen,  2) ; 
nan  *  ones  (indVarLen,  2); 
nan  *  ones  (indVarLen,  2); 
nan  *  ones  (indVarLen,  2); 

■  nan  *  ones  (indVarLen,  2) ; 
nan  *  ones  (indVarLen,  2); 
nan  *  ones  (indVarLen,  2); 
nan  *  ones  (indVarLen,  2); 
nan  *  ones  (indVarLen,  2) ; 
nan  *  ones  (indVarLen,  2); 
nan  +  ones  (indVarLen,  2) ; 
nan  *  ones  (indVarLen,  2) ; 
nan  *  ones  (indVarLen,  2); 

■  nan  *  ones  (indVarLen,  2); 

>  nan  *  ones  (indVarLen,  2); 

i  nan  *  ones  (indVarLen,  2); 

>  nan  *  ones  (indVarLen,  2); 

■  nan  *  ones  (indVarLen,  2); 

■  nan  *  ones  (indVarLen,  2); 


V,  INPUT  number  of  realizations  to  generate 
V.  number  of  realizations  accumulated 


see  later  calculation  of  pointing  error 


Initialize  graphics 


cmap  -  jet; 
invertBkgd  «  0 

sumBW  *=  0; 
figPat  *  1; 
figSum  -  2; 
cbarVert  -  0; 


V,  INPUT  color  map 

V.  INPUT  0:  figure  background  is  lowest  value  of  the  colormap; 

,,  j_.  «  highest 

INPUT  0:  summary  plot  in  color;  1:  in  black  and  white 
INPUT  figure  number  for  power  pattern 

and  summary  ,  ^  ... 

INPUT  0  for  horizontal  bar  below  pattern,  1  for  vertical  to  right 


V. 


V, 


V.  Set  window  positions 


0.15; 


cbarSize  • 
marginWid  -  8; 
jnarginHgt  -  44; 
screenSize  -  get 
screenWid  -  screenSize  (3) 
screenHgt  -  screenSize  (4); 
figPatWid  -  0.45  *  (screenSize 
figSumWid  -  screenSize  (3)  -  4 
if  cbarVert 
f igPatHgt 
figPatWid 
else 

f igPatHgt 
end 

figure  (figPat); 
figPat  -  gcf; 
set  (figPat 
1 posi tion 


V  colorbar  size  relative  to  pattern 
V,  width  margin  in  pixels 
%  height  margin  in  pixels 
(0,  ’screenSize*); 


(3)  -  4  *  marginWid); 

*  marginWid  -  figPatWid; 


figPatWid; 
figPatWid  * 


%  width  of  pattern  figure 
V.  width  of  summary  figure 

v.  height  of  pattern  figure 


(1  +  cbarSize) ; 


figPatWid  *  (1  +  cbarSize); 


create  and  position  pattern  figure 
update  in  case  figPat  couldn’t  be  created 


(marginWid 

screenHgt-marginHgt"f igPatHgt 
figPatWid 
f igPatHgt ) ,  ... 

•paperUnits ’ ,  ’inches') 
if  cbarVert 

set  (figPat,  ... 

’PaperOrientation’ ,  ’landscape’,  ... 

’ PaperPosition* ,  (0.5  0.5  10  7.5]); 

else 

set  (figPat,  — 
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•PaperOrientation' ,  ’portrait’ ,  ... 

* PaperPosi tion * ,  [0.5  0.5  7.5  10)); 

end 
elf ; 

figure  (figSum);  V.  create  and  position  summary  figure 

figSum  =  gcf;  V.  update  in  case  figSum  couldn't  be  created 

set  { f igSum,  . . . 

'position',  [screenWid-marginWid-figSumWid  ... 
marginHgt 
f igSumWid 

screenHgt-2*marginHgt] ) ; 

orient  tall; 
elf  ; 


Set  pattern  figure  colors 


if  invertBkgd 

bkgd  =  cmap  (size  (cmap,  1),  :);  4  take  background  from  high... 

else 

bkgd  -  cmap  (1,  :);  ?•  or  low  end  of  colormap 

end 

whitebg  (figPat,  bkgd);  «  set  pattern  colors 

set  (figPat,  ’color',  bkgd);  %  override  default  background  of  whitebg 


V, 


*  Set  summary  figure  colors  and  linestyles 


4 

whitebg  (figSum,  *w');  »  set  summary  colors 

set  (figSum,  ’color*,  ’w');  %  override  default  background  of  whitebg 

if  sumBW 

discrim  =  ’lineStyle*; 


discrimValue  =  *  —  *;  ‘4  line  styles 

discrimName  =  discrimValue;  %  names  for  legends  will  be  same  as  style  codes 

set  (figSum,  ...  ?.  black  color  forces  cycling  through  line  styles 

•DefaultAxesLineStyleOrder’ ,  discrimValue,  ... 

'DefaultAxesColorOrder' ,  [0  0  0) ) ; 


else 

discrim  -  'color*; 
set  (figSum,  ... 

•defaultAxesColorOrder*,  ’default’,  ... 
’defaultAxesLineStyleOrder’ ,  ’default’) ; 
discrimValue  -  get  (figSum,  ’defaultAxesColorOrder*); 
colorOrderHSV  »  rgb2hsv  (discrimValue) ; 
discrimValue  *»  num2cell  (discrimValue,  2); 
eolorNames  -  [ ' RYGCBMKW ' ) * ; 


discrimName  “  eolorNames  (1  +  ... 

mod  (round  (6  *  colorOrderHSV  (;,  1)),  6)); 
unsat  =  find  (colorOrderHSV  (:,  2)  <-  0.25); 
if  -isempty  (unsat) 

discrimName  (unsat)  «  eolorNames  (7  +  ... 
(colorOrderHSV  (unsat,  3)  >  0.5)); 

end 

discrimName  *  cellstr  (discrimName) ; 
clear  colorOrderHSV  eolorNames  unsat; 
end 


*4  reset  (if  previously  b&w,  for 
?.  example) 

9,  put  colors  in  array  of  RGB  coordinates 
V.  equivalent  HSV  coordinates 
V.  put  each  row  in  a  cell 
%  first  six  by  hue,  then  black  &  white 
if.  convert  H  to  an  integer  0..5 
‘4  then  index  into  eolorNames 
‘4  unsaturated  (gray)  colors 

%  threshold  V  to  0  (K)  or  1  (W) 

%  then  index  into  eolorNames 

%  put  each  letter  in  a  cell 


V.  Set  miscellaneous  common  properties 


set  ([figPat  figSum},  ... 

’invertHardCopy' ,  'off*,  ... 

’defaultTextFontSize' ,  8,  ... 

’defaultAxesFontSize' ,  8,  ... 

’toolbar’,  'none'); 
drawnow; 

clear  marginWid  marginHgt  screenSize  screenWid  screenHgt; 
clear  figPatWid  figSumWid  figPatHgt  bkgd; 

■4  Loop  over  independent  variable 
'4 

for  indx  «  1  :  indVarLen 

9.  Set  value  of  independent  variable 

V. 

eval  UindVarName  •  =  indVar  (indx);’]);  9.  set  variable  to  value 

•4  Code  may  be  inserted  below  to  transform  the  independent  variable  to 
V.  known  program  variables. 
stdTimeM  =  stdTimeMPS  *  le-3; 


Define  labels  for  substructures 

Here  we  construct  row  and  column  vectors  (corresponding  to  x  and  y, 
respectively)  that  tell  to  which  subaperture,  subarray,  or  element  a 


given  position  corresponds.  The  ; 
V.  numLY  =  2,  numMX  =  numMY  =  3,  and 

numLMX  =  numLX  *  numMX; 
numLMY  -  numLY  *  numMY ; 
numLMNX  -  numLMX  *  numNX; 
numLMNY  -  numLMY  *  numNY; 
nx  -  0  :  numLMNX  -  1; 

ny  =  (0  :  numLMNY  -  1) * ; 
mx  =  floor  (nx  /  numNX) ; 


;ample  vectors  are  for  numLX  = 
numNX  -  numNY  =  2. 

total  number  of  subarrays 

*4  total  number  of  elements 

•4  e.g.,  (0  1  2  3  4  5  6  7  8  9  10  11] 

■4  [0  1  2  3  4  5  6  7  8  9  10  11) 

(0011223344  5  5] 


?c'k 

VO'S 
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352  my  -  floor  <ny  /  numNY);  V,  [0  0  1  1  2  2  3  3  4455]’ 

353  lx  -  floor  (nx  /  (numMX  *  numNX));  V.  [0  0  0  0  0  0  1  1  1  1  1  1] 

354  ly  •=  floor  (ny  /  (numMY  *  numNY));  V.  [0  0  0  0  0  0  1  1  1  1  1  1]’ 

355 

356  V.  Calculate  ideal  (error-free)  excitation  magnitudes 

357 

358  excMagldl  «=  ones  (numLMNY,  numLMNX )  ;  V.  uniform  weighting 

359  if  diamond  V.  zero  every  other  element 

360  excMagldl  ■=  excMagldl  . 

361  .*  rem  (ones  (numLMNY,  1)  *  nx  +  ny  *  ones  (1,  numLMNX),  2); 

362  end 

363  excMagldl  -  excMagldl  /  sum  (excMagldl  (:).  A2)  ;  V.  normalize  to  unit  power 

364 

365  *?.  Calculate  ideal  (error-free)  excitation  phases 

366  V. 

367  V,  The  ideal  excitation  phases  are  those  that  produce  perfectly 

368  V.  constructive  interference  in  the  direction  of  the  steering  vector  at 

369  V.  the  reference  frequency.  This  implies  a  linear  phase  progression 

370  V, 

371  V,  excPhsIdl  -  -kO  (dx  sx  nx  +  dy  sy  ny)  +  const., 

372  V. 

373  V,  where  kO  (-  2  pi  fRef)  is  the  reference  wave  vector.  Quantization 

374  V.  prevents  the  array  from  achieving  this  flat  phase  front  for  all 

375  %  steering  angles.  However,  the  beamformer  simulated  here  mitigates 

376  l  effects  due  to  subaperture  quantization  by  adjusting  the  subarray 

377  %  delays,  which  is  effective  if  the  subarrays  are  quantized  at  a  finer 

378  V.  interval  than  the  subapertures.  Likewise,  it  compensates  for 

379  V.  subarray  quantization  with  the  element  phasers. 

380  V. 

381  V,  We  choose  the  constant  in  the  phase  progression  to  be  zero,  which 

382  V.  means  that  the  lowest  and  leftmost  components  (those  for  which  lx, 

383  v.  ly,  mx,  my,  nx,  or  ny  is  zero)  are  never  delayed,  while  the  highest 

384  Vi  and  rightmost  components  have  delays  that  depend  strongly  on  the 

385  V.  steering  vector. 

386  V, 

387  sx  -  sin  (steerPolar)  *  cos  (steerAzmth)  ;  V.  direction  cosines 

388  sy  *  sin  (steerPolar)  *  sin  (steerAzmth);  V.  for  steering 

389  sz  -  cos  (steerPolar);  *  used  much  later 

390  Time LX  -  -dx  *  sx  *  numNX  *  numMX  *  lx  /  c; 

391  TimeLY  -  -dy  *  sy  *  numNY  *  numMY  *  ly  /  c; 

392  TimeL  -  ones  (numLMNY,  1)  *  TimeLX  +  TimeLY  *  ones  (1,  numLMNX); 

393  if  qntTimeL  —  0  *  quantize? 

394  TimeL  -  round  (TimeL  /  qntTimeL)  *  qntTimeL;  V,  yes 

395  end; 

396  TimeMX  -  -dx  *  sx  *  numNX  *  mx  /  c; 

397  TimeMY  -  -dy  *  sy  *  numNY  *  my  /  c; 

398  TimeM  -  ones  (numLMNY,  1)  *  TimeMX  +  TimeMY  *  ones  (1,  numLMNX)  -  TimeL; 

399  if  qntTimeM  0 

400  TimeM  -  round  (TimeM  /  qntTimeM)  *  qntTimeM; 

401  end; 

402  TimeNX  -  -dx  *  sx  *  nx  /  c; 

403  TimeNY  -  -dy  *  sy  *  ny  /  c; 

404  PhseN  -  twopi  *  fRef  *  ... 

405  (ones  (numLMNY,  1)  *  TimeNX  +  TimeNY  *  ones  (1,  numLMNX)  -  TimeL  -  TimeM) ; 

406  if  qntPhseN  0 

407  PhseN  -  round  (PhseN  /  qntPhseN)  *  qntPhseN; 

408  end; 

409  excPhsIdl  *  twopi  *  fopr  *  (TimeL  +  TimeM)  +  PhseN; 

410  clear  TimeLX  TimeLY  TimeL  TimeMX  TimeMY  TimeM  TimeNX  TimeNY  TimeN  PhseN; 

411 

412  Vi  Prepare  transform  mapping 

413  Vi 

414  %  The  far-field  array  factor  is  the  Fourier  transform  of  the  complex 

415  Vi  excitations.  Considering  the  x  dimension  only  (the  operation  in  the 

416  4  y  dimension  is  analogous),  the  discrete  Fourier  transform  (DFT)  used 

417  4  later  calculates  the  far-field  array  factor  at  the  direction  cosine 

418  V,  cx  as 

419  4 

420  v.  tx-l 

421  V,  g  (cx  )  -  sum  exp  (-i  k  cx  dx  q)  e  , 

4  22  V.  q-0  q 

423 

424  v.  where  the  e(q]  are  the  complex  excitations  (zero-padded  if  tx 

425  V,  exceeds  the  array  size),  k  (-  2  pi  /  lambda)  is  the  operating  wave 

426  V.  vector,  and  lambda  (=  c  /  fOpr)  is  the  operating  wavelength.  The 

427  V.  argument  of  the  exponential  in  the  transform  may  be  written  -i  2  pi 

428  V.  (cx  /  lambda)  *  (dx  q) ,  where  cx  /  lambda  is  the  spatial  frequency 

429  v.  and  dx  q  is  the  spatial  coordinate.  The  direction  cosines  for  which 

4  30  V,  g  is  calculated  in  the  DFT  are 

431 

432  V.  lambda  p 

4  33  %  cx  - ,  P  -  0,  1,  .  .  .,  tx  -  1 

4  34  v.  P  dx  tx 

4  35 

436  v.  so  that  the  argument  of  the  exponential  is  -i  2  pi  p  q  /  tx.  We 

4  37  V.  have 

4  38 

439  V.  g  (cx  )  -  g  (cx  )  for  any  integer  p. 

440  v.  p+tx  p 

4  41  t 

442  •),  Given  g  (cx(p]),  the  array  factor  may  be  obtained  at  any  angle  using 

443 
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tx-l  /  p  dx  cx  \ 

g  (cx)  =  sum  S  I  —  - - ,  tx  I  g  (cx  ) 

p=0  \  tx  lambda  /  p 


A  where  the  geometric  progression 


1  tx-l 

S  (x,  tx)  =  —  sum  exp  (i  2  pi  q  x) 
tx  q=0 


1  exp  (i  2  pi  x  tx)  -  1 


tx  exp  (i  2  pi  x)  -  1 

is  an  interpolating  function,  but  this  formula  is  not  used  below. 
Incidentally,  note  that 

1  |  sin  (tx  pi  x)  | 

|S  (x,  tx)  |  =  —  I  -  I  / 

tx  |  sin  (pi  x)  | 


A  an  expression  that  often  appears  in  array  theory. 

A 

A  After  the  DFT  is  calculated,  the  program  maps  the  results  into  the 
A  region  of  cx-cy  space  where  the  direction  cosines  have  magnitude  1 
A  or  less,  tiling  as  necessary  to  fill  the  region.  The  portion  of 
A  that  region  for  which  cxA2  +  cyA2  <-  1  corresponds  to  visible  real 
A  space  (radiating  waves).  Later  processing  requires  a  border  of  at 
A  least  one  element  outside  the  visible  region.  This  section  prepares 
A  the  mapping. 


'A  values  of  p  (not  necessarily 
'A  integer)  for  which  cx  and  cy  are  1 
‘A  largest  integers  p  for  which 
A  cx,  cy  >  -1 

‘A  smallest  integers  p  for  which  cx, 

A  cy  >  (1  ♦  half  element  spacing) 


txLim  =  tx  *  dx  *  fOpr  /  c; 

tyLim  »  ty  *  dy  *  fOpr  /  c; 

txIndxLimMin  *  -floor  (txLim  )  -  1; 

tylndxLimMin  -  -floor  (tyLim  )  -  1; 

txIndxLimMax  -  floor  (txLim  +  0.5)  +  1; 

tylndxLimMax  -  floor  (tyLim  +  0.5)  +  1; 

•A  Extra  half  element  spacing  is  needed  only  for  surface  plots  with 
A  flat  shading;  see  graphics  code  below 
ax  -  txIndxLimMax  -  txIndxLimMin  +1;  A 

ay  “  tylndxLimMax  -  tylndxLimMin  +1;  A 

txlndx  -  txIndxLimMin  :  txIndxLimMax  ;  A 

tylndx  *  (tylndxLimMin  :  tylndxLimMax)';  A 

dirCosX  -  txlndx  /  txLim;  'A 

dirCosY  «=  tylndx  /  tyLim;  % 

dirCosShiftX  -  (txlndx  -  0.5)  /  txLim; 
dirCosShiftY  =  (tylndx  -  0.5)  /  tyLim;  % 

txlndx  -  txlndx  -  tx  *  floor  (txlndx  /  tx);  % 
tylndx  -  tylndx  -  ty  *  floor  (tylndx  /  ty) ;  A 

A  The  above  lines  accomplish  the  tiling  function  by  folding  the 
‘A  indices  into  the  interval  (0,  tx  -  1] 

tlndx  -  tylndx  *  ones  (1,  ax)  ...  ‘A  indices  to  elements  (one-based, 

+  ones  (ay,  1)  *  txlndx  *  ty  +  1;  A  column-ordered) 

clear  txIndxLimMin  tylndxLimMin  txIndxLimMax  tylndxLimMax  txlndx  tylndx; 


number  of  angle  samples  in  x 
and  y 

row  vector  of  indices 
column  vector 

corresponding  direction  cosines 
cx  and  cy 
A  shifted  direction  cosines  for  use 
with  flat  shading 
zero-based  indices  into  columns 
(x)  and  rows  (y)  of  DFT  results 


‘A  Prepare  far-field  angle  mapping 


V.  Physically,  the  array  factor  is  a  function  of  position  on  a 
V.  hemisphere.  The  direction  cosines  used  in  the  Fourier  transform  are 
V.  the  x  and  y  coordinates  of  points  on  the  unit  hemisphere.  Below  we 
•A  determine  the  region  of  the  transform  results  that  corresponds  to 
‘A  visible  space,  namely  cxA2  +  cyA2  <*  1,  and  calculate  the  z 
A  coordinates  of  points  in  visible  space  according  to 
A 

A  22  1/2 

A  cz  -  (1  -  cx  -  cy  )  ,  Re  cz  >-  0. 


A 

A  For  points  outside  visible  space,  cz  is  set  to  zero. 


A 

dirCosXMtx  -  ones  (ay,  1)  *  dirCosX; 
dirCosYMtx  =  dirCosY  *  ones  (1,  ax); 
radSqr  =  dirCosXMtx . A2  +  dirCosYMtx . A2 ; 
visBool  =  logical  (radSqr  <  1) ; 
dirCosZMtx  =  zeros  (ay,  ax); 

dirCosZMtx  (visBool)  =  sqrt  (1  -  radSqr  (visBool)); 


A  now  a  matrix 
A  squared  radius 

A  1  in  visible  space,  0  elsewhere 
'A  0  outside  visible  space 
A  positive  in  visible  space 


clear  radSqr; 


Loop  over  realizations 


for  rlzNura  =  1  :  numRlz 


",  Calculate  excitation  magnitudes  with  error 


excMagErr  =  1  +  stdAmplN 
err  =  1  +  stdAmplM 
excMagErr  =  excMagErr  . * 
err  =  1  +  stdAmplL 
excMagErr  =  excMagErr  .* 
excMag  «  excMagErr  .* 
clear  err  excMagErr; 


*  randn  (numLMNY,  numLMNX) ; 

*  randn  (numLMY  ,  numLMX  ) ; 
err  (my  +  1 ,  mx  +  1 ) ; 

*  randn  (numLY  ,  numLX  ) ; 
err  (ly  +  1,  lx  +  1)  ; 
excMagldl; 


A  element-level  error 
A  temporary  matrix 
A  subarray-level  error 
A  temporary 

A  subaperture-level  error 
A  actual  (with  error)  magnitude 


'A  Calculate  excitation  phases  with  error 


o*-  c 
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s, 

excPhsErr  *  stdPhseN  *  randn  (numLMNY,  numLMNX) ; 
err  -  stdTimeM  *  randn  (numLMY  ,  numLMX  ) 

*  twopi  *  fOpr; 

excPhsErr  =  excPhsErr  +  err  (my  +  1,  mx  +  1); 
err  =  (of sTimeL  +  stdTimeL  *  randn  (numLY, 

excPhsErr  =  excPhsErr  +  err  (ly  +1,  lx  +  1); 
excPhs  «=  excPhsErr  +  excPhsIdl; 
clear  err  excPhsErr; 


V.  element-level  error 
. .  ’■},  temporary  matrix 

‘A  (equivalent  phase) 

V.  subarray-level  error 
numLX) )  *  twopi  *  fOpr; 

V.  subaperture-level  error 
actual  (with  error)  phases 


V.  Assemble  complex  excitations 

exc  *=  excMag  .  *  exp  (-i  *  excPhs); 

V.  clear  excMag  excPhs; 

V.  Calculate  the  field  pattern  (array  factor) 

•*,  Here  the  DFT  is  calculated  and  the  result  rearranged  into  the 
V.  desired  region  of  direction  cosine  space.  The  element  factor  is 
V.  unity.  Mutual  coupling  is  ignored. 


g  -  f f t2  (exc,  ty,  tx); 
g  -  g  (tlndx) ; 

gSqr  -  real  (conj  (g)  .*  g); 
gMag  -  sqrt  (gSqr) ; 
clear  g; 


-A  first  element  is  zero  frequency 
?.  rearrange 
if,  squared  magnitude 
A  magnitude 


•f,  Determine  actual  beam  direction  by  fitting  the  transform 
‘l 

V.  This  first  of  two  methods  for  calculating  the  beam  pointing  vector 
V,  uses  the  information  in  the  Fourier  transform  of  the  excitations. 

V,  First  we  locate  the  element  of  g  with  the  largest  value.  (If  the 
V,  maximum  value  of  elements  is  obtained  by  more  than  one  element, 

V.  this  code  will  use  the  one  with  the  smallest  column  index  and  the 
V,  smallest  row  index  within  that  column.  Later  processing  will 
*  determine  whether  the  multiple  maxima  are  all  within  the  main 
V.  beam.)  To  estimate  the  location  of  the  maximum  of  the  underlying 
V.  continuous  function,  that  element  and  its  eight  neighbors  are 
V,  fitted  to  the  elliptic  paraboloid  [1] 
if, 

V.  1  2  12 

«*,  -Ux  +  Wxy  +  -Vy  +Xx  +  Yy  +  Z-  !g(x,  y)l 

'A  2  2 

if. 

V.  (in  the  direction  cosine  coordinate  system)  in  a  least-squares 
V.  sense.  We  employ  the  technique  of  QR  decomposition  to  find  the 
A  least-squares  solution.  Define  the  solution  vector 
V, 

V.  T 

a  -  [U  W  V  X  Y  Z] 

if, 

?.  of  length  M  -  6,  the  ordinate  (column)  vector  b  with  elements 

V. 

'I  b  -  I  g  (x  ,  y  )  I 

V.  i  i  i 

V, 

V,  of  length  N  “  9  (number  of  fitted  points),  and  the  design  matrix  A 
V,  having  rows 

•A  12  12 

V,  A  -I-x  xy  -y  x  y  1], 

V.  i,  :  2  i  i  i  2  i  i  i 

A 

%  where  the  (x(i],  y[i])  are  the  coordinates  of  the  points 

if,  neighboring  and  including  the  maximum  element.  QR  decomposition 

if.  of  A  factorizes  it  as 

V, 

V,  A  -  Q  R  , 

V.  where  Q  is  unitary  (O’  Q  =  eye)  and  R  is  upper  triangular.  (We 
V,  use  Matlab’s  economy  size  decomposition,  for  which  Q  is  M-by-N  and 
•!.  r  is  N-by-N.)  The  least-squares  solution  of 


A  a  *  b 


",  is  given  by 


a  -  R  Q*  b  . 

•j.  The  error  in  the  solution  depends  on  the  degree  to  which  the 
7,  values  of  Igl  depart  from  parabolic  form,  which  is  indicated  by 
the  reduced  chi-square 

2  12 
if.  chi  *  —  chi  , 

V,  nu  nu 

V.  where  nu  -  N  -  M  is  the  number  of  degrees  of  freedom  and 
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C  JFesf 


626  4  2 

627  4  chi  =  (A  a  -  b]  '  (A  a  -  b)  . 

628 

629  v.  The  covariance  matrix  for  the  solution  vector  is  normally  given  by 

630  4  the  matrix  inverse  of  the  curvature  matrix 

631 

632  4  alpha  =  A’  A  , 

633 

634  .  4  whose  elements  are  the  second  partial  derivatives  of  chi*2  with 

635  4  respect  to  the  elements  of  the  solution  vector  (2): 


636 

637  v.  -  1  d  d  2 

638  4  alpha  =  - - - chi 


639  *  mn  2d  a(m]  d  a(n] 

64  0 

641  V.  To  incorporate  the  degree  of  deviation  from  parabolic  form,  we 

64  2  V.  scale  the  covariance  matrix  by  the  reduced  chi-square,  as 

643  4 

644  V.  2  -12  -1 

64  5  4  C  -  chi  (A*  A)  -  chi  (R*  R)  . 

64  6  V.  poly  nu  nu 

647 

648  4  This  completes  the  least-squares  procedure. 

64  9  % 

650  4  Having  obtained  the  coefficients  of  the  best-fit  polynomial,  we 

651  4  locate  its  maximum.  The  coordinates  (plx,  ply)  of  the  maximum 

652  4  solve 

653  4 

654  4  /  u  w  \  /  plx  \  /  X  \  /  0  \ 

655  4  I  It  1+11*11; 

656  4  \  w  v  /  \  ply  /  \  Y  /  \  o  / 

657 

658  4  that  is, 

659  4 

660  4  /  plx  \  1/V-W\/X\ 

661  4  |  I  -  -  I  111* 

662  4  \  ply  /  d  \  -w  u  /  \  y  / 

663 

664  4  where 

665  4 

666  4  2 

667  4  D  -  W  -  U  V 

668 

669  4  is  the  discriminant  of  the  polynomial  (and  the  negative 

670  4  determinant  of  the  matrix).  If  the  maximum  so  found  lies  outside 

671  4  the  interpolation  region,  the  region  is  expanded  by  one  sample  in 

672  4  each  direction  and  the  least-squares  fit  is  repeated.  This  loop 

673  4  continues  until  the  interpolation  region  exceeds  a  certain  size  or 

674  4  a  satisfactory  maximum  is  found.  Assuming  a  maximum  has  been 

675  4  found,  the  covariance  matrix  for  plx  and  ply  is  calculated  next. 

676  4  We  first  form  the  derivative  matrix  or  Jacobian 

677  4 

678  4  d  (plx,  ply) 

679  4  J  - 

680  4  p  d  (a) 

681 

682  4  /  d  plx  d  plx  d  plx  d  plx  d  plx  d  plx  \ 

683 

684  4 

685 

686  4 

687 

688  4  XdUdwdVdXdYdZ/ 

689 

690  4  The  plx  derivatives  are 

691  4 

692  4  d  plx  V  X  -  w  Y 

693  4  - * - v 

694  4  d  u  DA2 

695 

696  4  d  plx  W-2  Y  +  U  V  Y  -  2  W  V  Y 

697  4  - - - 

698  4  d  w  D*2 

699 

700  4  d  plx  W  X  -  U  Y 

701  4  - - - w 

702  4  d  V  D*2 

703 

704  4  d  plx  V 

705  4  - -  - 

706  V.  d  X  D 

707 

708  d  plx  w 

709  4  - 

710  4  d  Y  D 

711  4 

712  4  and  the  ply  derivatives  follow  by  exchanging  U  for  V  and  X  for  Y 

713  4  everywhere.  The  covariance  matrix  for  plx  and  ply  is  simply 

714 

715  4  T 

716  4  C  =  J  C  J  - 


dU  dw  dV  dX  dY  dZ 
d  ply  d  ply  d  ply  d  ply  d  ply  d  ply 


717 


p  poly  P 


P 
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v. 

7,  If  no  satisfactory  maximum  was  found  by  the  above  procedure#  the 
V.  location  of  the  maximum  of  |g|  is  used.  A  covariance  matrix  is 
V,  fabricated  for  which  the  area  of  the  2  sigma  ellipse  equals  the 
v.  area  of  four  grid  squares  to  indicate  the  uncertainty  in  the 
V.  actual  location  of  the  maximum. 

•«.  (l)  d.  H.  von  Seggern,  _CRC  Standard  Curves  and  Surfaces_.  Boca 
V.  Raton,  FL:  CRC,  1993. 


7.  [2]  P.  R.  Bevington  and  D.  K.  Robinson,  _Data  Reduction  and  Error 
•i,  Analysis  for  the  Physical  Sciences_,  2nd  ed.  New  York,  NY: 

7,  McGraw-Hill,  1992,  pp.  121-125. 


gSqrMaxCol 

gSqrMaxCol 


[ gSqrMaxAct,  gSqrMaxRow]  =  max  (gSqr  .*  visBool); 
( gSqrMaxAct,  gSqrMaxCol]  =  max  (gSqrMaxAct); 
gSqrMaxRow  <=  gSqrMaxRow  (gSqrMaxCol); 
intRad  “  1; 
plOK  “  0; 

while  -plOK  &  intRad  <  4 

plNbrs  «  [-intRad  :  intRad}; 
plx  -  dirCosXMtx  (gSqrMaxRow  +  plNbrs, 
ply  -  dirCosYMtx  (gSqrMaxRow  +  plNbrs, 
plz  -  gMag  (gSqrMaxRow  +  plNbrs, 

plx  -  plx  <: ) ; 
ply  -  ply  (:); 
plz  -  plz  ( :) ; 

pDesMtx  -  [plx. *plx/2  plx. *ply  ply.*ply/2  ... 

plx  ply  ones (size (plx) )] ; 

[Q,  R]  -  qr  (pDesMtx,  0) ; 

pPoly  -  R  \  (O*  *  plz); 

dof  -  length  (plz)  -  length  (pPoly); 

chiSqr  *=  plz  -  pDesMtx  *  pPoly; 

chiSqr  -  chiSqr*  *  chiSqr; 

pPolyVar  «  inv  (R'  *  R)  *  chiSqr  /  dof; 

pVec  -  - [pPoly (1)  pPoly (2);  pPoly (2)  pPoly (3)) 

\  [pPoly (4) ;  pPoly (5) ] ; 
plOK  «... 

(pVec  (1)  >  dirCosX  (gSqrMaxCol  -  intRad)) 

&  (pVec  (1)  <  dirCosX  (gSqrMaxCol  +  intRad)) 

&  (pVec  (2)  >  dirCosY  (gSqrMaxRow  -  intRad)) 

&  (pVec  (2)  <  dirCosY  (gSqrMaxRow  +  intRad)) 

if  -plOK 

intRad  -  intRad  +  1; 
end 
end 

if  plOK 

plx  *  pVec  (1) ; 
ply  -  pVec  (2); 

pDet  -  pPoly  <2)A2  -  pPoly  (1) 
plxDer 


7,  maximum  values  and  their  rows 
V.  overall  maximum  and  column 
V,  corresponding  row 
V.  interpolation  radius 


plNbrs) ; 
plNbrs) ; 


4  neighbors 


gSqrMaxCol  +  plNbrs); 


4  design  matrix 

7,  now  pDesMtx  -  Q  *  R;  Q'  *  Q  ■ 
V.  solves  pDesMtx  *  pPoly  =  plz 
V.  degrees  of  freedom  in  the  fit 
%  deviations  only 
%  now  sum  of  squared  deviations 
4  pDesMtx’  *  pDesMtx  -  R’  *  R 
4  find  critical  point 

4  is  interpolated  point 
l  inside  neighborhood? 

%  (pathological  cases  can 
7,  place  it  outside) 


%  keep  interpolated  point 


%  determinant 


pPoly  (3); 

[[I  -1] * (pPoly ( [3  2] ) . *pPoly ( (4  5] ) ) ‘pPoly (3) /pDet 
[1  1  -2] * (pPoly ( [2  1  2) ) ,*pPoly( [2  3  3] ) . *pPoly ( (5  5  4]))/pDet 
pPoly ( [4  51 ) ) *pPoly (2) /pDet 


U  -1] * (pPoly ( [2  1]) 
pPoly (3) 

-pPoly (2) 

0  ]  '  /  pDet; 

plyDer  -  [[1  -1] * (pPoly ( [2  3]) 

(1  1  -2]MpPoly((2  3 
[1  -1] * (pPoly([l  21) 
-pPoly (2) 
pPoly (1) 

0  )*  /  pDet; 

[plxDer;  plyDer]  *  pPolyVar 


4  derivatives 

*pPoly ( [5  4] ) ) *pPoly (2) /pDet 
2 ] ) . *pPoly ( ( 2  1  1} ) .*pPoly( [4  4  5]))/pDet 
‘pPoly ( [5  4 ) ) ) *pPoly ( 1 ) /pDet 


plVar 
else 

plx  *  dirCosXMtx  (gSqrMaxRow,  gSqrMaxCol); 
ply  -  dirCosYMtx  (gSqrMaxRow,  gSqrMaxCol); 
plVar  -  diag  (1  ./  (pi  *  [txLim  tyLim].A2)); 
end 

clear  intRad  plOK  plNbrs  plz  pDesMtx  Q  R  pPoly; 
clear  dof  chiSqr  pPolyVar  pVec  pDet  plxDer  plyDer; 


[plxDer;  plyDer]*;  4  covariance  matrix 

4  interpolated  point  is  outside 
4  use  location  of 
4  actual  maximum 
4  2  sigma  area  «  4  grid  squares 


7.  Determine  actual  beam  direction  by  fitting  the  excitation  phases 

%  The  second  method  for  calculating  the  pointing  vector  uses  the 
%  excitation  magnitudes  and  phases,  not  the  transform.  For  brevity, 
7,  define 

Delta  (x,  y)  =  k  n  dx  x  +  k  m  dy  y  +  theta 
'}.  mn  mn 

7,  where  the  theta  (m,n)  are  the  excitation  phases,  and  let  e[m,nl 
denote  the  excitation  magnitudes.  When  we  express  the  power 
pattern  as 


2 

lg  (x,  y) I  =  sum  sum  sum  sum  e  e 

ml  nl  m2  n2  ml,nl  m2,n2 


ml ,  nl 


)) 

m2,n2 


eye 


exp  [i  (Delta 


Delta 
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2 

=  sum  e  +  2  sum  sum'  e  e 

m,n  m,  n  ml,nl  m2,n2  ml,nl  m2,n2 

*  cos  (Delta  -  Delta  )  , 

ml , nl  m2 , n2 


where  the  primed  sum  is  over  distinct  pairs  (ml,  nl)  and  (m2,  n2), 
we  see  that  the  maximum  occurs  where  the  cosine  contributions  are 
largest.  If  the  phase  front  is  nearly  flat,  the  arguments  of  the 
cosines  will  be  small  for  (x,  y)  near  the  direction  of  phase  front 
propagation.  To  fourth  order  in  the  arguments. 


Ig  (x,  y) I 


(sum  e  ) 
m,  n  mn 


2  sum  sum  e 

ml,nl  m2,n2  ml,nl 


m2 ,  n2 


/  1  2 

|  -  (Delta  -  Delta  ) 

\  2  ml,nl  m2,n2 

1  4  \ 

+  —  (Delta  -  Delta  )  I  - 

24  ml,nl  m2,n2  / 


Keeping  terms  to  only  second  order,  we  are  motivated  to  find  the  x 
and  y  (implicit  in  Delta [mn])  that  minimize 


chi  =  sum  sum  e  e  (Delta  -  Delta  )  , 

2  ml,nl  m2,n2  ml,nl  m2,n2  ml,nl  m2,n2 

where  the  subscript  2  denotes  the  second-order  truncation.  Taking 
derivatives  with  respect  to  x  and  y  and  rearranging  yields  the 
normal  equations 

/  0  \  /  nl-n2  \ 

|  |  «=  sum  sum  eel  I 

\  0  /  ml,nl  m2,n2  ml,nl  m2,n2  \  ml-m2  / 

/  \ 

*  |  (nl-n2  ml-m2}  Pi  +  theta  -  theta  I  , 

\  ml,nl  m2,n2  / 

where  Pi  -  k  [dx*x  dy*y] ’ .  Here  the  sums  contain  a  total  of 
(numLMNXA2  numLMNYA2)  terms,  which  may  be  of  the  order  of  one 
million  for  a  typical  array.  To  reduce  this  number,  we  transform 
the  least-squares  problem  to  an  equivalent  but  simpler  problem. 
First,  the  above  normal  equations  may  be  rewritten  with  the  column 
vector  [nl-n2  ml -m2] »  replaced  by  [nl  ml]*,  which  may  be  seen  by 
separating  the  column  vector  into  two  sums  and  exchanging  (ml,nl) 
with  (m2,n2)  in  one  of  the  sums.  Second,  the  row  vector  may  be 
separated  into  two  sums  to  obtain 


/  0  \ 

|  |  =  (  sum  e 

\  0  /  m,  n  i 


/  n  \ 

)  sum  e  I  |  (  [n  m]  Pi  +  theta 

nm,nm,n\m/  m 


V,  /  /  n  \  \ 

•i  —  |  sum  e  I  II  sum  e  (  (n  mj  Pi  +  theta  ) 

V,  \  m,n  mn  \  m  /  /  m,n  mn  mn 

V. 

V,  Based  on  the  second  term,  we  define 
‘i 

$  sum  e  (  [n  m]  Pi  +  theta  ) 

$  m,  n  mn  mn 

V.  Delta  *  - -  - 

sum  e 

«j>  m,  n  mn 

V.  Finally,  we  may  rewrite  the  normal  equations  as 


/  0  \ 

|  0  |  —  sum 
\  0  /  m,  n 


/ 

e  I 
m,n  \ 


\ 

1  (  (n  i 

/ 


1)  Gamma  +  theta  ) 


where  Gamma  =  [k*dx*x  k*dy*y  Delta]’  and  the  third  row  follows 
from  the  definition  of  Delta.  These  normal  equations  contain  only 
(numLMNX  numLMNY)  terms,  nominally  on  the  order  of  1000.  They 
find  the  plane  k  n  dx  x  +  k  m  dy  y  +  Delta  that  best  fits 
-theta (m,n)  in  a  weighted  least-squares  sense.  The  solution  is 
found  using  OR  decomposition  of  the  design  matrix,  which  has  rows 
(n  m  1 ] . 


\\  The  error  in  the  solution  is  determined  not  by  the  deviations  of 
%  the  -theta [mn]  from  the  best- fit  plane  nor  by  the  deviations 
Delta  [ml , nl ]  -  Delta  [m2, n2]  appearing  in  chiA2  earlier,  for  the 
V,  solution  is  exactly  the  power  pattern  maximum  to  second  order  in 
the  cosine  arguments.  However,  the  error  does  depend  on  the 
V,  fourth  and  higher  powers  of  the  cosine  arguments.  So  motivated, 
"}.  we  consider  the  fourth-order  merit  function 
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2 

V,  chi  «*  sum  sum  e  e  (Delta 

V,  4  ml,nl  m2,n2  ml,nl  m2,n2  ml,nl 

/I  2  \ 

♦II---  (Delta  -  Delta  )  I 

\  12  ml,nl  m2 , n2  / 

and  observe  that  chi4A2  »<  chi2A2  everywhere.  We  interpret  the 
V.  difference  chi2A2  -  chi4A2  as  indicative  of  the  error  in  the 
V.  solution,  and  we  scale  the  covariance  matrix  by  that  amount.  The 
covariance  matrix  used  is  the  inverse  of  the  curvature  matrix  for 
V.  chi2A2;  that  curvature  matrix  is 

V.  2 

",  alpha  =  k  sum  sum  e  e 

tf,  ml,nl  m2,n2  ml,nl  m2,n2 

V,  /  (nl-n2 ) dx  \ 

V,  *  |  |  [  (nl-n2)dx  (ml-m2)dy  ]  . 

V.  V  (ml-m2)dy  / 

V. 

V,  Because  we  construct  the  design  matrix  for  the  simpler 
V,  least-squares  problem,  we  must  construct  alpha  explicitly. 

V.  However,  this  can  be  accomplished  by  analytically  expanding  the 
V,  differences  and  factoring  the  sums. 

Vi 

phsX  -  ones  (numLMNY,  1)  *  nx  *  dx  *  (twopi  *  fppr  /  c) ; 
phsY  -  ny  *  ones  (1,  numLMNX)  *  dy  *  (twopi  *  fOpr  /  c) ; 
phsX  -  phsX  (:) ; 
phsY  «  phsY  ( : )  ; 
excMagV  -  excMag  (:); 

desMtx  -  (phsX  phsY  ones (numLMNX* numLMNY, 1) 1  ... 

.*  (sqrt  (excMagV)  *  ones  (1,  3)); 

IQ,  R]  -  qr  (desMtx,  0);  V  now  desMtx  «  Q  *  R  and  Q'  *  Q  »  eye  [conj. 
excPhsWgt  -  -excPhs  (:)  .*  sqrt  (excMagV); 

P2Vec  -  R  \  (Q*  *  excPhsWgt) ;  *  desMtx  *  p2Vec  -  excPhsWgt 

p2x  -  p2Vec  (1); 
p2y  -  p2Vec  (2); 

DeltaPhs  -  (phsX  phsY  ones (numLMNX * numLMNY, 1) )  *  p2Vec  +  excPhs  (;); 

sumExcMagDeltal  -  excMagV  .*  DeltaPhs; 

sumExcMagDelta2  *  sumExcMagDeltal  .*  DeltaPhs; 

sumExcMagDelta3  -  sumExcMagDelta2  .*  DeltaPhs; 

sumExcMagDelta4  «  sumExcMagDelta3  .  *  DeltaPhs; 

sumExcMagDeltaO  «  sum  (excMagV); 

sumExcMagDeltal  -  sum  (sumExcMagDeltal); 

sumExcMagDelta2  -  sum  <sumExcMagDelta2) ; 

sumExcMagDelta3  -  sum  (sumExcMagDelta3) ; 

sumExcMagDelta4  -  sum  (suraExcMagDelta4 ) ; 

chiSqrRed  -  (  2  *  sumExcMagDelta4  *  sumExcMagDeltaO  ... 

-  8  *  sumExcMagDelta3  *  sumExcMagDeltal  . . . 

+  €  *  sumExcMagDelta2  *  sum£xcMagDelta2  )  /  12; 

chiSqrRed  -  max  (0,  chiSqrRed);  V,  in  case  of  roundoff  error 

excMagPhs  -  excMagV*  *  (phsX  phsY] ; 
crvMtx  -  2  *  sumExcMagDeltaO  *  (phsX  phsY]’  ... 

*  ( (phsX  phsY]  .♦  (excMagV  excMagV])  ... 

-  2  *  excMagPhs  *  *  excMagPhs; 
p2Var  -=  chiSqrRed  *  inv  (crvMtx); 

clear  phsX  phsY  excMagV  desMtx  Q  R  excPhsWgt  p2Vec  DeltaPhs; 

clear  sumExcMagDeltaO  sumExcMagDeltal  sumExcMagDelta2  sumExcMagDelta3 

clear  chiSqrRed  excMagPhs  crvMtx; 

V,  Construct  pointing  vector 
V. 

V.  Above  we  constructed  two  pointing  vectors  by  different  methods. 

V.  The  method  of  fitting  the  transform  is  robust  even  for  large 
V.  errors  but  limited  by  the  transform  resolution.  On  the  other 
V.  hand,  the  method  of  fitting  the  excitation  phases  is  independent 
V,  of  transform  resolution  but  accurate  only  for  small  errors, 

V,  approaching  the  exact  solution  as  the  phase  errors  decrease.  We 
V,  wish  to  obtain  a  single  pointing  vector  for  subsequent  use,  and 
V,  for  this  purpose  we  form  a  weighted  average.  Specifically,  each 
•«.  vector  is  weighted  by  the  inverse  of  the  area  of  its  covariance 
■*,  ellipse,  which  is  pi  times  the  determinant  of  the  covariance 
v.  matrix.  Similarly,  a  single  covariance  matrix  is  obtained  by 
V,  weighing  each  covariance  matrix  by  the  square  of  the  pointing 
",  vector  weights,  normalized  to  avoid  effectively  halving  the 
v.  covariance  matrix  when  the  two  incoming  matrices  are  nearly  equal. 

areal  =  det  (plVar) ; 

area2  =  det  (p2Var) ; 

areaTot  =  areal  +  area2; 

wghtl  =  area2  /  areaTot; 

wght2  -  areal  /  areaTot; 

px  =  wghtl  *  plx  +  wght2  *  p2x; 

py  -=  wghtl  *  ply  +  wght2  *  p2y; 

pVar  =  (wghtlA2  *  plVar  +  wght2A2  *  p2Var)  /  (wght2A2  +  wghtlA2) ; 
clear  areal  area2  areaTot  wghtl  wght2; 
clear  plx  ply  plVar; 
clear  p2x  p2y  p2Var; 

V.  Calculate  peak  power  density  and  pointing  error 


2 

-  Delta  ) 

m2 ,  n2 


transpose] 


sumExcMagDelta4 ; 
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ft 

V.  To  avoid  inaccuracies  due  to  interpolation,  the  peak  power 
V.  density  is  obtained  by  explicitly  evaluating  the  Fourier 
V.  transform  at  the  pointing  vector.  The  pointing  error  is 
V.  straightforwardly  calculated  from  the  dot  product  of  the  steering 
V,  and  pointing  vectors.  However,  an  optional  second  method  is 
V.  coded  that  makes  use  of  the  pointing  vector  covariance  matrix  to 
V.  calculate  the  uncertainty  (standard  deviation)  of  the  pointing 
‘v.  error  due  to  uncertainty  in  the  pointing  vector.  If  this 

uncertainty  is  desired,  also  uncomment  lines  elsewhere  that  refer 
V,  to  errPointUnc  and  errPointUncS . 

V. 


l); 

exp  (-i  *  twopi  *  (fOpr  /  c)  ... 

*  dx  *  ones  (numLMNY,  1)  +  nx  ... 

*  dy  *  ny  *  ones  (1,  numLMNX) ) )  . *  exc; 
sum  (gSqrMax  (:)); 

real  (conj  (gSqrMax)  *  gSqrMax) ; 

(1  -  pxy2); 


pxy2  =  px/'2  +  py*2 
peakVisb  =  (pxy2 
if  peakVisb 
gSqrMax  = 

*  (  px 

+  py 
gSqrMax  - 
gSqrMax  = 
pz  -  sqrt 
if  0 

pCrs  -  [  0  -sz  sy 
sz  0  -sx 

-sy  sx  0  J  *  (px  py  pz)*; 
pCrsDer  -  (  -sy*px/pz  -sz-sy*py/pz 
s z+sx  *px/pz  sx*py/pz 
-sy  sx  ] ; 

pCrsDer  *  pVar  *  pCrsDer'; 
sqrt  (pCrs*  *  pCrs); 


complex  field 


phaser  sum 
power 


ft  INPUT  0  to  skip  std  dev, 
ft  cross  product 


%  rows  correspond  to 
ft  components  of  pCrs; 
%  columns,  to  pVec 
%  covariance  matrix 
ft  magnitude 


pCrs*  /  pCrsMag; 
pCrsMagDer  *  pCrsVar 

trace  (pCrsVar)  /  3; 


V.  derivative  exists 


derivative  doesn't  exist 
average  of  principal  variances 


pCrsVar 
pCrsMag 
if  pCrsMag  > 
pCrsMagDer 

pCrsMagVar  =  pCrsMagDer  *  pCrsVar  *  pCrsMagDer*; 
else 

pCrsMagVar 
end 

errPoint  *=  asin  (pCrsMag)  ; 

errPointUnc  *  abs  (1  /  sqrt  (1  -  pCrsMagA2))  *  sqrt  (pCrsMagVar); 
clear  pCrs  pCrsDer  pCrsVar  pCrsMag  pCrsMagDer  pCrsMagVar; 
else 

errPoint  -  acos  (min  (1,  ...  ft  dot  product  for  error;  min 

sx  *  px  +  sy  *  py  +  sz  *  pz));  ft  prevents  roundoff  problems 

errPointUnc  -  nan; 
end 

ejse  ft  maximum  is  invisible 

gSqrMax  =  gSqrMaxAct ; 
px  <*  nan; 
py  =  nan; 
pz  -  nan; 
errPoint  -  nan; 
errPointUnc  -  nan; 
end 

beamPowerDB  -  10  *  loglO  (gSqrMax); 
clear  pxy2; 


ft  Determine  main  beam  region 


V,  The  angular  domain  of  the  main  beam  is  constructed  starting  with 
V.  the  maximum  element.  The  largest  neighboring  element  is  added  on, 
V,  followed  by  the  largest  neighbor  of  either  point,  and  so  on.  This 
V,  accretion  continues  until  any  neighbor  of  the  largest  element  on 
ft  the  main  beam  border  exceeds  the  element  added  previously. 

V,  Effectively,  elements  are  added  with  values  descending  from  the 
%  peak  until  an  opportunity  to  ascend  is  reached.  All  visible 
ft  elements  outside  of  the  main  beam  are  declared  to  be  in  the 
ft  sidelobes . 


beamWidLvl  “  gSqrMax  /  2; 
adjc  «  {-ay-1  -ay  -ay+1  -1 

V.  Note:  We  must  include  the  diagonal 
ft  correctly  descend  a  structure  such  as 
beamBool  =  logical  (zeros  (ay,  ax)}; 
brdrLen  *=  1; 
brdrlndx  =  (gSqrMaxCol  -  1)  *  ay  +  gSqrMaxRow; 
brdrVal  =  0; 

beamBool  (brdrlndx)  =  1; 

adjclndx  =  brdrlndx  +  adjc'; 

adjclndx  =  adjclndx  (visBool  (adjclndx)); 

adjcVal  =  gSqr  (adjclndx) ; 

beamDepth  =  inf; 

beamVisb  =  1; 

capVisb  =  ( gSqrMaxAct  >  beamWidLvl); 

while  max  (adjcVal)  <=  beamDepth 
brdrLen  =  brdrLen  -  1; 
brdrlndx  =  brdrlndx  (1  ;  brdrLen); 
brdrVal  =  brdrVal  (1  :  brdrLen); 
beamBool  (adjclndx)  -  ones  (size  (adjclndx)); 
for  adjcPtr  =  1  :  length  (adjclndx) 

pos  -  sum  (brdrVal  <-  adjcVal  (adjcPtr)); 
brdrlndx  =  (brdrlndx (1 :pos)  adjclndx (adjcPtr) 
brdrVal  =  {brdrVal ( 1 :pos)  adjcVal (adjcPtr) 
brdrLen  =  brdrLen  +  1; 


0.5  0.9). 

V.  build  main  beam  in  Boolean  variable 

!(.  start  with  maximum 
any  value  will  do  here 
main  beam  begins  with  maximum 
7,  and  neighbors 

that  are  visible 
get  values 

get  the  loop  started 

7.  usually  true  unless  resolution  is  too  low 
are  the  new  neighbors  all  downhill? 
yes;  remove  element  from  border 


■i,  add  neighbors  to  main  beam 
V.  and  to  border 
ft  ordered  least  to  greatest 
brdrlndx (pos+1 tbrdrLen) ) ; 
brdrVal (pos+1 :brdrLen)  ); 


ft  power  level  where  beam  width  is  measured 
1  ay-1  ay  ay+1);  ft  relative  indices  of  neighbors 

neighbors  in  order  to 
[1  0.4 
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'«*,  pick  largest  element  from  border 
'■},  neighbors  of  chosen  element 
V.  eliminate  invisible  points 
V.  were  some  invisible? 

V.  yes;  clear  flag 

are  we  below  the  threshold? 

1  no;  the  cap  is  partially  invisible 


V.  use  only  new  elements 
V.  and  get  their  values 

V,  closed  contour  at  beamWidLvl? 
sidelobe  region 

?,  duplicate  maximum  outside  beam? 
%  no;  the  beam  is  identified 


end 

beamDepth  *=  brdrVal  (brdrLen)  ; 
adjclndx  «  brdrlndx  (brdrLen)  +  adjc; 
adjclndx  «  adjclndx  (visBool  (adjclndx)); 
if  length  (adjclndx)  <  length  (adjc) 
beamVisb  =  0; 

if  beamDepth  >=  beamWidLvl 
capVisb  =0; 
end 
end 

adjclndx  «=  adjclndx  (-beamBool  (adjclndx)); 
adjcVal  =  gSqr  (adjclndx); 
end 

capClosed  -  capVisb  &  (beamDepth  <  beamWidLvl); 
sideBool  »  visBool  &  -beamBool; 
beamlndx  -  find  (beamBool); 

if  max  (max  (gSqr  (sideBool)))  <  gSqrMaxAct 
beamExist  -  1; 
if  -peakVisb 

disp  (’Warning:  The  main  beam  peak  is  invisible;  some  calculations' ) ; 
disp  ('  may  return  NaN.'); 
elseif  -capVisb 

disp  (’Warning:  The  beam  width  contour  of  the  main  beam  is  partially*}; 
disp  {’  invisible;  some  calculations  may  return  NaN.*); 
elseif  -capClosed 

disp  ('Warning:  The  main  beam  is  insufficiently  deep  for  obtaining*); 
disp  ('  its  width;  some  calculations  may  return  NaN.’); 
elseif  -beamVisb 

disp  ('Warning:  The  main  beam  is  partially  invisible;  some*); 
disp  (’  calculations  may  return  NaN.*); 
end 
else 

beamExist  «  0;  'A  yes;  the  beam  is  ambiguous 

disp  ('Warning:  The  main  beam  is  not  identifiable;  some  *); 
disp  (*  calculations  will  return  NaN.*); 
px  -  nan; 
py  -  nan; 
pz  *  nan; 
errPoint  *  nan; 
beamVisb  -  0; 
beamlndx  -  U; 
peakVisb  -0; 
capVisb  -  0; 
capClosed  *»  0; 
beamDepth  *  nan; 
sideBool  -  visBool; 
end 

if  beamDepth  --  0 
beamDepthDB  "  -inf; 
else 

beamDepthDB  -  10  *  loglO  (beamDepth  /  gSqrMax) ; 
end 

clear  adjc  beamBool  brdrLen  brdrlndx  brdrVal  adjclndx  adjcVal  beamDepth  adjcPtr  pos; 


%  strike  earlier  results 


treat  visible  space  as  sidelobes 


%  avoid  warning  message 


V,  Determine  main  beam  width  and  roll 

V, 

■A  The  analysis  of  the  beam’s  width  and  roll  is  conducted  using  a 
if,  stereographic  projection,  for  which  projections  of  great  circles 
V,  intersect  at  the  same  angles  as  the  great  circles  on  a  sphere. 

V,  (See  the  comments  in  the  plotting  section  below  for  details.) 

V.  This  property  allows  us  to  obtain,  in  the  limit  of  a  narrow  beam, 

V,  the  correct  roll  angle  and  the  beam  widths  along  two  orthogonal 
%  great  circles. 

?. 

!?,  The  actual  calculations  are  based  on  fitting  the  half-power  contour 
V.  of  the  main  beam  to  an  ellipse.  First  the  contour  is  obtained  in 
V,  direction  cosine  space,  then  the  coordinates  are  transformed  to 
V.  stereographic  coordinates.  The  contour  is  fitted  to  the  conic 
V,  section 


•f,  12  12 

V,  -  U  x  +  Wxy+-Vy  +  Xx  +  Yy+Z*=0 


'■!.  using  a  simple  algorithm  that  minimizes  the  algebraic  distance  as 
V,  follows.  Define  the  design  matrix  D  to  have  rows 

l  12  12 

v.  D  «I(-x)(xy)(-y)x  y  1), 

V,  i,  :  2  i  i  i  2  i  i  i 

v.  where  the  (x(i),  y  [  i ) )  are  the  points  along  the  contour,  and  let  the 

coefficient  vector  be 


V.  T 

a  *=  [U  w  V  X  Y  Z]  . 

V,  The  algebraic  distance  between  a  point  and  a  conic  section  is  the 
V.  left-hand  side  of  the  conic  section  equation,  so  that  the  distance 
V.  between  a  point  i  along  the  contour  and  the  ellipse  described  by  the 
V,  vector  a  is  simply  D(i,:]a.  We  seek  the  minimum  of  the  sum  of 
V.  squared  algebraic  distances,  which  is  just  ||D  a||A2,  subject  to  the 
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7  constraint  ||a||*2  -  1.  We  therefore  introduce  the  constrained 
7,  objective  function 

2 

V,  E  -  |  |D  all  -  lambda  (Halt  -  1) 

T  T  T 

7,  =  a  D  D  a  -  lambda  (a  a  -  1) 

where  lambda  is  a  Lagrange  multiplier.  The  minimum  is  found 
7.  analytically  to  occur  when 

V,  T 

7,  D  D  a  =  lambda  a  , 

7.  which  is  an  eigenvalue  equation.  The  desired  coefficient  vector,  a, 
V.  corresponds  to  the  minimum  eigenvalue. 

I  Using  the  coefficients  of  the  best-fit  ellipse,  we  now  calculate 
7,  the  beam  characteristics.  First,  a  sign  change  is  applied  to  the 
7  coefficients  if  necessary  to  force  U  (and  therefore  V)  to  be 
V.  negative.  For  convenience,  we  rewrite  the  conic  section  as 
7 

V,  IT  T 

7  -pAp  +  Bp  +  Z«0, 

7,  2 

7 

7  where  p  =  [x;  y] , 

7  /  U  W  \ 

7  A  -  I  I  / 

7  \  W  V  / 

7 

7  and  B  «  [X;  Y] .  We  first  find  the  center  of  the  ellipse  in  order 
7  to  draw  it  later.  Replacing  p  with  p  +  pi  in  the  conic  section 
7  yields 
7 

7  IT  T  T 

7  -  pi  A  pi  +  (p  A  +  B  )  pi  +  Z1  =  0  , 

7  2 


where 

IT  T 

Z1  *  -  p  Ap  +  B  p+Z 
2 


7  is  defined  for  later  use.  When  p  coincides  with  the  center,  the 
7  linear  term  vanishes;  therefore,  p  solves 


7  A  p  +  b  =  0  . 

7 

7  We  next  find  the  roll  angle,  which  is  conceptually  defined  as 
7  follows,  using  spherical,  not  stereographic,  coordinates.  If  the 
7  ellipse  center  is  not  at  boresight,  rotate  it  (and  the  antenna 
7  pattern)  to  boresight  along  the  great  circle  connecting  the  two. 

7  The  angle  from  the  great  circle  with  azimuth  0  to  the  great  circle 
7  along  the  beam’s  major  axis  (direction  of  maximum  width)  is  the 
7  roll  angle.  Alternatively,  construct  the  great  circle  connecting 
7  the  ellipse  center  and  boresight.  The  roll  angle  is  the  sum  of 
7  two  angles,  the  angle  from  the  great  circle  with  azimuth  0  to  the 
7  constructed  great  circle  and  the  angle  from  the  constructed  great 
7  circle  to  the  great  circle  along  the  beam's  major  axis.  Now  the 
7  roll  angle  so  defined  is  merely  the  apparent  orientation  of  the 
7  major  axis  when  viewed  in  the  stereographic  projection.  In  a 
7  (stereographic)  coordinate  system  rotated  by  that  angle,  the 
7  off-diagonal  element  of  A  (the  coefficient  W)  vanishes;  therefore, 
7  we  seek  the  coordinate  system  that  diagonalizes  A.  Replacing  p 
7  with  R  p2  in  the  original  conic  section  gives 


ITT  T 

-  p2  (R  A  R)  p2  +  (B  R)  p2  +  Z  =  0  . 
2 


V,  The  new  quadratic  coefficient  RAT  A  R  will  be  diagonal  if  the 
7  columns  of  R  are  the  eigenvectors  of  A.  The  new  diagonal  elements 
7  U2  and  V2  become  the  eigenvalues,  which  are  explicitly 


U  +  V  /  2  /  U  -  V  \2  \l/2 

U2  - - +  |  W  +  |  - I  I  and 

2  \  \  2  /  / 


U  +  V 

V2  - - 

2 


/  2  /  U  -  V  \2  U/2 

j  W  +  |  - i  I 

\  \  2  /  / 


7  The  eigenvalue  with  the  smaller  magnitude  (U2  above,  since  U  and  V 
7  are  negative)  corresponds  to  the  major  axis.  Therefore,  the 
7  corresponding  eigenvector  points  along  the  direction  of  the  major 
7  axis;  the  other  eigenvector,  along  the  minor  axis.  The  roll  angle 
7  is  obtained  from  the  two  components  of  the  major  axis;  explicitly, 
7  it  satisfies 
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7  2  W 

7  tan  (2  roll)  - - . 

U  -  V 

"  We  also  use  the  eigenvectors  to  draw  the  ellipse  later.  Last,  the 
V,  eigenvalues  yield  the  major  and  minor  full  widths  at  half  maximum 
7.  of  the  main  beam  as 

7.  /  Z1  \l/2  /  Z1  \l/2 

|-8  —  |  and  !  -8  —  |  , 

\  U2  /  \  V2  / 

7,  respectively.  As  these  were  derived  in  the  stereographic 

■?.  projection,  a  factor  of  (1  +  cos  polar)  is  applied  to  obtain  the 

7,  approximate  widths  in  real  angles. 


if  capClosed 

V. 

7  Construct  the  contour 


caplndx  -  beamlndx  (gSqr  (beamlndx)  >-  beamWidLvl);  V.  elements  at  or  above  level 
capBool  -  logical  (zeros  (ay,  ax)); 
capBool  (caplndx)  «  ones  (size  (caplndx)); 
adjc  -  {ay  ay+1  1  -ay+1  -ay  -ay-1  -1  ay-1);  7 

dirlndx  -  1;  * 

intlndx  -  caplndx  (1);  ‘4 

dirlndxSt  -  0;  7 


intlndxSt  »  0 
capContX  =  [] 
capContY  *  () 
while  capBool 


clockwise  in  matrix  row-column  coordinates 
initial  index  into  adjc 
initial  index  of  interpolation  center 
get  loop  started 


(intlndx  +  adjc  (dirlndx)) 


intlndx  +  adjc  (dirlndx); 


7.  empty  contour  coordinates 


7.  next  element  is  inside  cap? 

V,  keep  moving  until  edge  is  reached 


%  back  at  starting  point? 

7  no;  get  magnitude 
7  looking  across  row  or  column? 
7  across  column? 

7  yes;  interpolate  in  x 


intlndx 
end 

while  (intlndx  —  intlndxSt)  I  (dirlndx  —  dirlndxSt) 
adjclnc  -  abs  (adjc  (dirlndx)); 
if  (adjclnc  —  1)  I  (adjclnc  «*=  ay) 
if  adjclnc  —  ay 

capContYl  -  dirCosYMtx  (intlndx); 

capContX 1  -  dirCosXMtx  (intlndx)  +  (beamWidLvl  -  gSqr  (intlndx))  ... 

*  (dirCosXMtx  (intlndx  +  adjc  (dirlndx))  -  dirCosXMtx  (intlndx))  ... 

/  (gSqr  (intlndx  +  adjc  (dirlndx))  -  gSqr  (intlndx)); 

else  7  across  row 

capContXl  -  dirCosXMtx  (intlndx);  7  interpolate  in  y 

capContYl  -  dirCosYMtx  (intlndx)  +  (beamWidLvl  -  gSqr  (intlndx))  ... 

*  (dirCosYMtx  (intlndx  +  adjc  (dirlndx))  -  dirCosYMtx  (intlndx))  - 


/  (gSqr  (intlndx  +  adjc  (dirlndx))  -  gSqr 

end 

if  isempty  (capContX)  ! 

capContX  -  capContXl;  ! 

capContY  *  capContYl; 
intlndxSt  =  intlndx; 
dirlndxSt  -  dirlndx; 

elseif  (capContX  (length  (capContX))  —  capContXl) 

&  (capContY  (length  (capContY))  capContYl) 
capContX  -  [capContX;  capContXl); 
capContY  =  (capContY;  capContYl); 

end 


(intlndx) ) ; 

7  first  point? 

7  yes;  store  it 

7  remember  starting  point 

—  7  duplicate?  (e.g., 

7  element  equals  level) 
7  no;  append  it 


end 

dirlndx  -  dirlndx  +  1; 
if  dirlndx  >  length  (adjc) 
dirlndx  ■  1; 
end 

adjclnc  -  abs  (adjc  (dirlndx)); 
if  capBool  (intlndx  +  adjc  (dirlndx) ) 
intlndx  -  intlndx  +  adjc  (dirlndx) ; 
dirlndx  -  dirlndx  -  length  (adjc)  /  2  +  1; 
if  (adjclnc  —  1)  I  (adjclnc  --  ay) 
dirlndx  «  dirlndx  +  1; 
end 

if  dirlndx  <  1 

dirlndx  =  dirlndx  +  length  (adjc); 
end 
end 

end  V,  while 


7  (no  diagonal  interpolation) 
7  next  direction 
7  cycle 


7  no;  get  magnitude  of  step 
7  next  element  is  inside  cap? 

7  yes;  becomes  new  interpolation  center 
7  reverse,  then  ahead  one  increment 
7  stepped  in  row  or  column? 

7  yes;  ahead  an  extra  increment 
7  (useless  to  look  back  diagonally) 

7  cycle 


contour  complete 


7  Fit  an  ellipse  in  stereographic  coordinates 

capContZ  -  sqrt  (1  -  capContX. A2  -  capContY. *2 ) ; 
capContXS  =  capContX  ./  (1  +  capContZ); 
capContYS  =  capContY  ./  (1  +  capContZ); 

cDesMtx  =  (capContXS . *capContXS/2  capContXS . *capContYS 
capContYS . *capContYS/2  capContXS  ... 
capContYS  ones (size (capContXS) ) ] ; 
cDesMtx  -  cDesMtx*  *  cDesMtx; 

capContX  -  capContX  ([ 1 : length (capContX)  1]); 
capContY  =  capContY  ([ 1 : length (capContY)  1]); 
capContZ  -  capContZ  ( [1 : length (capContZ)  1)); 

(cEigVec,  cEigVal)  -  eig  (cDesMtx); 

(cEigValMin,  cEigValMinldx)  -  min  (diag  (cEigVal)); 
cPoly  *=  cEigVec  (:,  cEigValMinldx) ; 
cPoly  -  -cPoly  *  sign  (cPoly  < 1 ) ) ; 


7  to  stereographic  coords 

7  to  form  design  matrix 
7  for  least-squares  fit 

7  done 

7  close  contour  for  plotting 


7  eigenvectors  and  -values 
7  minimum  eigenvalue 
7  and  matching  vector 
7  to  have  negative  eigenvalues  below 
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v. 

V,  Obtain  beam  characteristics  from  ellipse  coefficients 


cHessMtx  =  lcPoly(l)  cPoly(2) 

cPoly (2 J  cPoly  <3)1; 
cDervMtx  =  l cPoly (4) 

cPoly (5) ] ; 

capCenter  =  -cHessMtx  \  cDervMtx; 

capConst  =  cPoly  (6)  +  cDervMtx’  *  capCenter  /  2; 

(cEigVec,  cEigVal ]  =  eig  (cHessMtx) ; 

[cEigVal,  cEigValOrd)  =  sort  (diag  (cEigVal) ) ; 
cEigVec  -  cEigVec  (:,  cEigValOrd); 
if  cEigVec  (1,  2)  ==  0 
roll  =  pi  /  2; 
else 

roll  =  atan  (cEigVec  (2,  2)  /  cEigVec  (1,  2)); 
end 

roll  -  roll  -  pi/2  +  azmthOffst; 

roll  -  roll  -  ceil  (roll  /  pi)  *  pi  +  pi/2  -  azmthOffst; 
hpbw  «  sqrt  (-8  *  capConst  .  /  cEigVal); 


%  Hessian  matrix; 

V.  second  derivatives 
V.  first  derivative  matrix 

Vi  ellipse  center 

Vi  new  constant  coefficient 

V,  ascending  order 
V.  corresponding  order 
V.  special  case? 

V.  angle  to  major  axis 
V.  make  -pi/2  < 

V.  roll  +  azmthOffst  <-  pi/2 
V,  two-vector 


V,  Construct  the  fitted  ellipse  for  plotting  (in  stereographic  coordinates) 
% 

theta  -  (0  :  10  /  360  :  1)  *  twopi; 
cFit  -  capCenter  *  ones  (size  (theta)) 

+  0.5  *  (cEigVec  *  diag  (hpbw)  ... 

*  l sin (theta);  cos (theta) ]) ; 
cFitR2  -  sum  (cFit.A2) 


calculate  some 
points  along  the 
fitted  ellipse 

squared  radius  in  stereographic  coords 

z  direction  cosine 

undo  stereographic  projection 


V,  undo  widths,  too 


cFitZ  -  (1  -  cFitR2)  ./  (1  +  cFitR2); 
cFitX  -  cFit  (1,  :)  .*  d  +  cFitZ); 
cFitY  *  cFit  (2,  :)  .*  (1  +  cFitZ); 
hpbw  «  hpbw  *  (1  +  pz); 
hpbwMjr  -  hpbw  (2); 
hpbwMnr  «  hpbw  (1); 
else 

roll  -  nan; 
hpbwMjr  -  nan; 
hpbwMnr  -  nan; 
end 

clear  capBool  adjc  dirlndx  intlndx  dirlndxSt  intlndxSt  adjclnc  capContXl  capContYl; 
clear  capContXS  capContYS  cDesMtx  cEigVec  cEigVal  cEigValMin  cEigValMinldx  cPoly; 
clear  cHessMtx  cDervMtx  capCenter  capConst  cEigValOrd  hpbw; 
clear  theta  cFit  cFitR2; 


V»  Calculate  power  in  visible  space,  main  beam,  and  sidelobes;  main 
V>  beam  and  sidelobe  solid  angles;  and  average  sidelobe  level 


%  These  calculations  involve  integrals  over  the  hemisphere  or  portions 
%  of  it.  The  integrals  are  carried  out  in  direction  cosine 
%  coordinates  by  multiplying  the  integrand  by  the  appropriate 
Vi  Jacobian. 


%  The  integrals  are  evaluated  using  the  midpoint  approximation,  for 
V.  which  the  starting  point  is  the  Taylor  series  expansion  of  g  (x,  y) 
V,  to  second  order: 

V, 

V.  g  (a  +  u,  b  +  v) 


V,  «lg  +  ug+vg+-ug  +uvg  +-vg  I 

\  x  y  2  xx  xy  2  yy  /a,b 

V, 

V.  Then 


Vi  /  c/2  /  d/2 

V.  !  du  |  dv  g  (a  +  u,  b  +  v)  - 

/-C/2  /-d/2 

l  /  1  2  2  \ 

V,  =  |  c  d  g  +  —  cd(cg+dg)l 

V,  \  24  xx  yy  /a,b 

V, 

V  The  midpoint  approximation  keeps  the  first  term  and  neglects  the 
quadratic  terms.  The  maximum  amount  neglected  is 

c  d  /  2  |  I  2  1  l  \ 

V,  - |c  Ig  (a,  b)  1  +  d  Ig  (a,  b)||, 

24  \  l  xx  1  l  yy  I  / 


\«.  which  we  use  as  the  error  estimate  for  each  interior  point  of  the 
V.  integration,  approximating  the  second  derivatives  with  scaled  second 
differences. 


V.  To  the  interior  error  is  added  an  estimate  of  the  error  due  to 
V.  finite  sampling  at  the  integral  limits;  the  estimate  is  half  the 
V.  value  of  the  integrand  at  the  outermost  samples. 

V.  (The  error  estimate  calculations  have  been  commented  out  for 
V.  speed.) 


V.  visbEdgelndx  *  visBool; 

V.  VisbEdgelndx  (2  :  ay  -  1,  2  :  ax  -  1) 
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1453  V,  <  visbEdgelndx  (2  :  ay  -  1,  3  :  ax  )  S  visbEdgelndx  (1  :  ay  -  2,  2  :  ax  -  1)  ... 

1454  7,  4  visbEdgelndx  (2  :  ay  -  1,  1  :  ax  -  2)  &  visbEdgelndx  (3  :  ay  ,2  :  ax  -  1)  ); 

1455  V,  visbEdgelndx  *  find  (visBool  -  visbEdgelndx); 

1456  V.  beamEdgelndx  ■  zeros  {ay,  ax); 

1457  V.  beamEdgelndx  (beamlndx)  =  ones  (size  (beamlndx)  )  ; 

14  58  ‘.l.  sideEdgelndx  ■  beamEdgelndx; 

1459  V.  beamEdgelndx  (2  ;  ay  ►  1,  2  :  ax  -  1)  **  . . . 

1460  t  (  beamEdgelndx  (2  :  ay  -  1,  3  :  ax  )  &  beamEdgelndx  {1  :  ay  -  2,  2  :  ax  -  1 )  ... 

1461  &  beamEdgelndx  (2  ;  ay  -  1,  1  :  ax  -  2)  4  beamEdgelndx  (3  :  ay  ,  2  :  ax  -  1)  ); 

14  62  V.  beamEdgelndx  (beamlndx)  ■=  1  -  beamEdgelndx  (beamlndx); 

1463  V.  beamEdgelndx  =  find  (beamEdgelndx); 

14  64  V.  sideEdgelndx  (2  :  ay  -  1,  2  ;  ax  -  1)  =  — 

14  65  Vi  (  sideEdgelndx  (2  :  ay  -  1,  3  :  ax  )  I  sideEdgelndx  ( 1  :  ay  -  2,  2  :  ax  -  1 )  ... 

14  66  'I,  I  sideEdgelndx  (2  :  ay  -  1,  1  :  ax  -  2)  I  sideEdgelndx  (3  :  ay  ,  2  :  ax  -  1)  ); 

14  67  V.  sideEdgelndx  (beamlndx)  =  zeros  (size  (beamlndx)); 

14  68  V.  sideEdgelndx  «  find  (sideEdgelndx); 

1469  areaFact  -  zeros  (ay,  ax);  V.  areal  factor  for  integrating: 

1470  areaFact  (visBool)  «...  V.  Jacobian  (1  /  cos  polar) 

1471  1  ./  (txLim  *  tyLim  *  dirCosZMtx  (visBool));  V.  and  grid  spacing 

1472  areaFactLim  -  2  /  sqrt  (txLim  *  tyLim);  l  edge  points  may  exceed  this 

1473  tooBig  -  find  (areaFact  >  areaFactLim);  V.  arbitrary  limit;  force  those 

1474  areaFact  (tooBig)  -  ones  (size  (tooBig))  *  areaFactLim;  «  that  do  to  comply 

1475  V  sldAng  -  sum  (areaFact  (:)); 

1476  '4  sldAngErrRel  «  abs  (sldAng  /  (2  *  pi)  -  1); 

1477  %  areaFactUnc  ■  zeros  (ay,  ax); 

1478  %  areaFactUnc  (2  :  ay  -  1,  2  :  ax  -  1)  -  . . .  t  error  estimate  based  on  second- 

1479  %  (  abs  (diff  (areaFact  (2  :  ay  •  1,  2)')  ...  %  order  Taylor  expansion 

1480  7,  +  abs  (diff  (areaFact  ( 2  :  ax  -  1)  ,  2)  )  )  /  24; 

1481  V.  areaFactUnc  ([1  ay),  :)  -  areaFactUnc  ([2  ay-1),  :>;  %  assume  errors  at  outer  edge  equal 

1482  V,  areaFactUnc  {:,  (1  ax])  -  areaFactUnc  (:,  [2  ax-1]);  V.  those  of  nearest  neighbors 

1483  intgrnd  *=  areaFact  .*  gSqr;  V.  to  integrate  gSqr 

1484  V,  intgrndUnc  -  zeros  (ay,  ax);  l  error  estimate  for  integrand 

1485  V.  intgrndUnc  (2  :  ay  -  1,  2  :  ax  -  1)  »  — 

1486  l  (  abs  (diff  (intgrnd  (2  :  ay  -  1,  :)\  2)')  ... 

1487  V.  +  abs  (diff  (intgrnd  ( : ,  2  :  ax  -  1)  ,  2)  )  )  /  24; 

1488  7,  intgrndUnc  ([1  ay],  :)  -  intgrndUnc  ((2  ay-1],  :);  V.  assume  errors  at  outer  edge  equal 

1489  Vi  intgrndUnc  (:,  [1  ax))  -  intgrndUnc  (:,  {2  ax-1]);  ?.  those  of  nearest  neighbors 

1490  powerVisb  -  sum  (intgrnd  (:)); 

1491  Vi  powerVisbUnc  -  sum  (intgrndUnc  (:)); 

1492  directivityDB  -  10  *  loglO  (4  +  pi  *  gSqrMax  /  powerVisb); 

1493  if  beamExist 

1494  7,  sldAngMain  -  sum  (areaFact  (beamlndx));  V,  main  beam  solid  angle 

1495  7.  sldAngMainUnc  -  sum  (areaFactUnc  (beamlndx))  +  sum  (areaFact  (beamEdgelndx))  /  2; 

1496  powerMain  -  sum  (intgrnd  (beamlndx));  V.  power  in  the  main  beam 

1497  V.  powerMainUnc  »  sum  (intgrndUnc  (beamlndx))  +  sum  (intgrnd  (beamEdgelndx))  /  2; 

1498  powerSide  -  sum  (intgrnd  (sideBool));  V  power  in  the  sidelobes 

1499  V  powerSideUnc  -  sum  (intgrndUnc  (:))  -  sum  (intgrndUnc  (beamlndx)) — 

1500  V,  +  sum  (intgrnd  (sideEdgelndx))  /  2; 

1501  sldAngSide  -  sum  (areaFact  (sideBool));  V.  sidelobe  equivalent  solid  angle 

1502  V.  sldAngSideUnc  -  sum  (areaFactUnc  (sideBool))  -  sum  (areaFactUnc  (beamlndx))  ... 

1503  V  +  sum  (areaFact  (sideEdgelndx))  /  2; 

1504  else 

1505  V,  sldAngMain  -  nan; 

1506  V.  SldAngMainUnc  »  nan; 

1507  powerMain  -  nan; 

1508  V.  powerMainUnc  -  nan; 

1509  powerSide  -  nan; 

1510  V.  powerSideUnc  -  nan; 

1511  sldAngSide  =  nan; 

1512  V,  sldAngSideUnc  «  nan; 

1513  end 

1514  powerMainVisbDB  -  10  *  loglO  (powerMain  /  powerVisb);  V  ratio  of  main  beam  to  visible  power 

1515  **  powerMainVisbUnc  -  powerMainUnc  /  powerVisb  ... 

1516  V.  +  powerMain  *  powerVisbUnc  /  powerVisb*2; 

1517  powerVisbSideDB  -  10  *  loglO  (powerVisb  /  powerSide);  V.  ratio  of  visible  to  sidelobe  power 

1518  V.  powe rVi sbS ideUnc  -  powerVisbUnc  /  powerSide  ... 

1519  V.  +  powerVisb  *  powerSideUnc  /  powerSide^; 

1520  powerMainSideDB  -  10  *  loglO  (powerMain  /  powerSide);  V.  ratio  of  main  beam  to  sidelobe  power 

1521  V.  powe rMainS ideUnc  -  powerMainUnc  /  powerSide  ... 

1522  V.  +  powerMain  *  powerSideUnc  /  powerSide'' 2; 

1523  powe r S ideAvgDB  *  10  *  loglO  (powerSide  /  (sldAngSide  *  gSqrMax));  V.  average  sidelobe  power 

1524  V.  powerSideAvgUnc  =  powerSideUnc  /  sldAngSide  ... 

1525  +  powerSide  *  sldAngSideUnc  /  sldAngSide'' 2; 

1526  clear  intgrnd  areaFact  areaFactLim  tooBig; 

1527  V.  clear  intgrndUnc  areaFactUnc; 

1528 

1529  V.  Locate  nearest  and  largest  sidelobes 

1530 

1531  V.  We  wish  to  identify  the  sidelobe  closest  in  angle  to  the  main  beam 

1532  V.  and  the  sidelobe  with  the  largest  peak  power.  We  first  find  the 

1533  "  local  maxima  in  the  sidelobe  region,  then  for  each  we  determine 

1534  V.  the  possible  ranges  for  actual  distance  from  the  main  beam  and 

1535  V.  peak  power.  (The  uncertainties  arise  from  the  discrete  sampling 

1536  V.  of  the  array  factor.)  Using  the  ranges  we  select  those  peaks  that 

1537  V.  could  possibly  be  the  closest  or  largest.  For  each  of  these 

1538  candidates  a  more  precise  location  and  peak  power  is  computed  by 

1539  V.  interpolating  over  neighboring  data.  Finally,  based  on  these 

1540  V.  results,  the  closest  and  nearest  sidelobes  are  identified. 

1541 

1542  if  peakVisb 
154  3 

1544  V,  Find  local  maxima  (sidelobe  peaks),  using  discrete  differences 
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1554 

1555 

1556 

1557 

1558 

1559 

1560 

1561 

1562 

1563 

1564 

1565 

1566 

1567 

1568 

1569 

1570 

1571 

1572 

1573 

1574 

1575 

1576 

1577 

1578 

1579 

1580 

1581 

1582 

1583 

1584 

1585 

1586 

1587 

1588 

1589 

1590 

1591 

1592 

1593 

1594 

1595 

1596 

1597 

1598 

1599 

1600 
1601 
1602 

1603 

1604 

1605 

1606 

1607 

1608 

1609 

1610 
1611 
1612 

1613 

1614 

1615 

1616 

1617 

1618 

1619 

1620 
1621 
1622 

1623 

1624 

1625 

1626 

1627 

1628 

1629 

1630 

1631 

1632 

1633 

1634 

1635 

1636 


V.  to  approximate  derivatives.  The  differences  are  formed  from  the 
V.  magnitude  of  the  array  factor,  not  the  squared  magnitude;  since 
7.  the  behavior  should  already  be  parabolic  near  peaks,  squaring 
V.  would  produce  fourth-order  behavior  and  make  second-order 
V,  interpolation  less  accurate. 


V.  Key  to  the  variables  below: 

7.  x  first  differences  in  x 

V.  Y  first  differences  in  y 

XX  second  differences  in  x 

7.  YY  second  differences  in  y 

7,  X2  first  differences  in  x  with  double  step 

7,  XY  cross  differences  in  x  and  y 

V.  XC  nonzero  where  first  difference  in  x  changes  sign 

■?.  YC  nonzero  where  first  difference  in  y  changes  sign 


gMagX  =  gMag  (:,  2: ax)  -  gMag  <:,  l:ax-l); 
gMagY  -  gMag  {2: ay,  :)  -  gMag  <l:ay-l,  :); 
gMagXX  =  {zeros {ay, 1)  (gMagX  (:,  2:ax-l) 

gMagYY  -  [ zeros (1 , ax) ;  (gMagY  <2:ay-l,  :) 

gMagX 2  -  [zeros (ay, 1)  (gMag  (:,  3:ax  ) 

gMagX Y  -  [ zeros  (1, ax) ;  (gMagX2(3:ay  ,  :) 

gMagXC  -  [zeros (ay, 1)  (gMagX  (:,  2:ax-l)  . 
gMagYC  -  [zeros (1, ax) ;  (gMagY  (2:ay-l,  :)  . 
sllndx  -  find  {  sideBool  — 
t  gMagXC  t  gMagYC 

&  (gMagXX  <  0)  &  (gMagYY  <  0)  ... 

&  (gMagXY. A2  -  gMagXX  .*  gMagYY  <  0)  ); 
if  isempty  (sllndx) 
sllndx  =  [] ; 
slNrstDist  =  nan; 


-  gMagX  (:,  l:ax-2))  zeros { ay, 1 }] ; 

-  gMagY  (l:ay-2,  :));  zeros (1, ax) ] ; 

-  gMag  (:,  l:ax-2))/2  zeros (ay, 1) ] ; 

-  gMagX2(l:ay-2,  :))/2;  zeros ( 1, ax) ] ; 

*  gMagX  (:,  l:ax-2)  <  0)  zeros (ay, 1) ) ; 

*  gMagY  (l;ay-2,  :)  <  0);  zeros (1, ax) ] ; 

%  identify  sidelobe  points  where 

7,  first  derivatives  change  sign, 

V;  second  derivatives  are  negative,  an 
%  discriminant  is  negative 
V.  none  found? 


slNrstPowrDB  -=  nan; 
slNrstVec  =  nan  ♦  ones  (1,  3); 
slLgstDist  *  nan; 
slLgstPowrDB  ■  nan; 
slLgstVec  -  nan  *  ones  (1,  3); 

se  %  found  some  peaks 


7.  Compute  possible  ranges  of  distances  and  powers;  identify 
7,  candidates  for  closest  and  largest  peaks 


toward  *  sign  (px  -  dirCosXMtx  (sllndx))  *  ay  ... 

+  sign  (py  -  dirCosYMtx  (sllndx)); 
slCosDistMax  -  ... 

px  *  dirCosXMtx  (sllndx  +  toward)  . . . 

+  py  *  dirCosYMtx  (sllndx  +  toward)  . . . 

+  pz  *  dirCosZMtx  (sllndx  +  toward); 

slCosDistMin  «  ... 

px  *  dirCosXMtx  (sllndx  -  toward)  ... 

+  py  *  dirCosYMtx  (sllndx  -  toward)  ... 

+  pz  *  dirCosZMtx  (sllndx  -  toward); 

slNrstBool  -  (slCosDistMax  >=  max  (slCosDistMin) ) ; 
slPowr  -  (gMag  (sllndx)  ... 

+  (  -gMagXX  (sllndx)  ... 

+  2  *  abs  (gMagXY  (sllndx))  ... 

-  gMagYY  (sllndx)  )  /  8).A2; 
slLgstBool  -  (slPowr  >=  max  (gSqr  (sllndx))); 
slCandlndx  -  sllndx  (slNrstBool  I  slLgstBool); 
numCand  *  length  (slCandlndx) ; 

7. 

7.  Interpolate  powers  and  locations  for  candidates 


V,  index  increment  to  neighbor 
7.  closer  to  pointing  vector 
7.  cosine  of  maximum  possible 
7.  angle  between  pointing 
7.  vector  and  each  peak; 

7,  dot  product 
7  likewise  for  minimum 
%  possible  angle 


7.  true  if  peak  might  be  the  closest 
7.  estimate  largest  possible 
7.  interpolated  power  by  adding 
7.  an  error  estimate  based  on  the 

7,  differences  computed  above 

7.  true  if  peak  might  be  the  largest 
v.  candidates  for  closest  and  largest 
7.  number  of  candidates 


six  -  zeros  (numCand,  1) ; 
sly  ~  zeros  (numCand,  1); 
slPowr  =*=  zeros  (numCand,  1)  ; 
for  slPtr  *  1  :  numCand 

slCol  *  ceil  (slCandlndx  (slPtr)  /  ay) ; 

si  Row  -=  slCandlndx  (slPtr)  -  (slCol  -  1)  *  ay; 

intRad  =1; 


^  allocate  space  for  x  and 

%  y  direction  cosines 

V.  and  powers 

V.  loop  through  candidates 

V,  separate  index  into  column 

7.  and  row  indices 

7.  interpolation  radius 


slOK  =  0; 
while  -slOK 
slNbrs 
slxFit 
slyFit 
slzFit 
slxFit 
slyFit 
slzFit 
slDesMtx 


&  intRad  <  3 


[ -intRad  ;  intRad) ; 

dirCosXMtx  (slRow  +  slNbrs,  slCol  +  slNbrs); 

dirCosYMtx  (slRow  +  slNbrs,  slCol  +  slNbrs); 

gMag  (slRow  +  slNbrs,  slCol  +  slNbrs) ; 

slxFit  { ; ) ; 
slyFit  (:); 
slzFit  (:); 

-  (slxFit. ‘slxFit/2  slxFit. ‘slyFit  slyFit. ‘slyFit/2 
slxFit  slyFit  ones (size (slxFit) )) ; 

[0,  R]  =  qr  (slDesMtx,  0) ; 
slFoly  -  R  \  (O'  *  slzFit) ; 

slVec  -  - [slPoly (1)  slPoly<2);  slPoly(2)  slPoly<3)] 

\  (slPoly (4) ;  slPoly(5)]; 


7i  initialize 

7.  no  answer  yet  but  too  early  to  bail? 
7.  offsets  to  neighbors 


x,  y,  and  z  coordinates  of 
neighbors  (using  magnitude, 
not  power,  for  z) 


design  matrix 


now  slDesMtx  =  Q  *  R;  Qf  *  0  =  eye 
solves  slDesMtx  *  slPoly  =  slzFit 
7.  find  critical  point 


slOK  =  ... 
(slVec 
&  (slVec 
&  (slVec 
&  (slVec 
if  -slOK 
intRad  - 


(1)  >  dirCosX  (slCol  -  intRad))  . 

(1)  <  dirCosX  (slCol  +  intRad))  . 

(2)  >  dirCosY  (slRow  -  intRad) )  . 

(2)  <  dirCosY  (slRow  +  intRad) ) ; 

intRad  +1; 


is  interpolated  point 
inside  neighborhood? 
(pathological  cases  can 
place  it  outside) 

outside  neighborhood? 
yes;  cast  a  wider  net 


end 
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end  interpolation  attempts 
interpolation  successful? 
yes;  keep  interpolated  point 


end 

if  slOK 

six  (slPtr)  =  slVec  < 1 > ; 
sly  (slPtr)  =  slVec  (2) ; 

slPowr  (slPtr)  -  ...  interpolate  power 

((six (slPtr) . *slx(slPtr) /2  six (slPtr) . *sly (slPtr)  ... 

sly (slPtr) . *sly (slPtr) /2  slx(slPtr)  sly(slPtr)  1J  *  slPoly).A2; 


else 

six  (slPtr)  *  dirCosXMtx  (slCandlndx  (slPtr)); 
sly  (slPtr)  =  di rCosYMtx  (slCandlndx  (slPtr) ); 
slPowr  (slPtr)  =  gSqr  (slCandlndx  (slPtr)); 
end 
end 

Select  closest  and  largest  peaks 


interpolation  failed 
V.  use  grid  location  of 
V.  sampled  maximum 
",  use  sampled  power 

V.  end  of  loop  through  candidates 


slz  -  sqrt  (1  -  six . A2  -  sly . *2) ; 
slDist  -  acos  (px  *  six  +  py  *  sly  +  pz  *  slz); 

IslNrstDist,  slNrstlndx]  -min  (slDist); 

slNrstPowrDB  -  10  *  loglO  (slPowr  (slNrstlndx)  /  gSqrMax) ; 
slNrstVec  -  ( six (slNrstlndx)  sly (slNrstlndx)  slz (slNrstlndx) ) ; 
[slLgstPowrDB,  slLgstlndx)  -max  (slPowr); 
slLgstPowrDB  -  10  *  loglO  (slLgstPowrDB  /  gSqrMax); 
slLgstDist  -  slDist  (slLgstlndx) ; 

slLgstVec  -  (six (slLgstlndx)  sly (slLgstlndx)  slz (slLgstlndx) ) ; 
end 

else  4  peak  is  invisible 

sllndx  -  {); 
slNrstDist  -  nan; 

SlNrstPowrDB  -  nan; 
slNrstVec  -  nan  *  ones  (1,  3); 
slLgstDist  -  nan; 
slLgstPowrDB  *  nan; 
slLgstVec  -  nan  *  ones  (1,  3); 
end 

clear  gMagX  gMagY  gHagXX  gMagYY  gMagX2  gMagXY  gMagXC  gMagYC; 

clear  toward  slCosDistMax  slCosDistMin; 

clear  slNrstBool  slLgstBool  slCandlndx  numCand; 

clear  six  sly  slz  slPowr  slPtr  slCol  slRow  intRad  slOK; 

clear  slNbrs  slxFit  slyFit  slzFit  slDesMtx  Q  R  slPoly  slVec; 

clear  slDist  slnrstindx  slLgstlndx; 


V.  z  direction  cosines 
V.  angular  distances 

V.  smallest  distance 
V.  and  corresponding  power 
V.  keep  vector  for  plotting 
4  largest  power 
V,  converted  to  dB 
V,  and  corresponding  distance 
4  keep  vector 


In  order  to  calculate  running  means  and  standard  deviations  of  n 
realizations,  we  accumulate  the  mean  and  variance 


4  Record  characteristics 

4 

V. 

4 
4 
4 
4 
4 
4 
4 
4 
4 


M  -  -  sum  x 
n  n  m-1  m 


V  - - sum  (x 

n  n  -  1  m=l  m 


4  using  the  updating  formulas 


x  -  M 


n  -  2  1  2 

v  - - V  +  -  (x  -  M  )  . 

n  n  -  1  n-1  n  n  n-1 

The  running  mean  is  simply  M(n] ,  and  the  standard  deviation  is 
sqrt  (V{n)).  M  is  accumulated  in  the  first  column  of  a  matrix;  V 
in  the  second.  This  method  is  more  immune  to  roundoff  error  than 
accumulating  the  sums  of  values  and  squares  (1). 

4  If  the  beam  is  invisible,  the  excitations  and  array  factor  are 
7,  saved  to  an  automatically-named  file. 

7,  (i)  n.  J.  Higham,  ^Accuracy  and  Stability  of  Numerical  Algorithms_ 
4  Philadelphia,  PA;  SIAM,  1996,  pp.  12-13. 

if  beamVisb  I  (numRlz  —  1} 
if  isnan  (numAcc  (indx) ) 


1716 

1717 

numAcc  (indx)  «  1; 
pxS 

( indx. 

1) 

px; 

1718 

pxS 

(indx. 

2) 

= 

0; 

1719 

pys 

(indx, 

1) 

- 

py; 

1720 

pyS 

(indx. 

2) 

= 

0; 

1721 

pzS 

(indx. 

1) 

* 

pz; 

1722 

pzS 

(indx. 

2) 

■ 

0; 

1723 

errPointS 

(indx. 

1) 

- 

errPoint ; 

1724 

errPointS 

(indx. 

2) 

- 

0; 

1725 

4  errPointUncS 

(indx. 

1) 

- 

errPointUnc, 

1726 

4  errPointUncS 

(indx. 

2) 

- 

0; 

1727 

beamPowerDBS 

(indx. 

1) 

- 

beamPowerDB, 

1728 

beamPowerDBS 

(indx. 

2) 

0; 
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1729 

beamDepthDBS 

(indx,  1) 

- 

beamDepthDB; 

1730 

beamDepthDBS 

(indx,  2) 

= 

0; 

1731 

hpbwMj  rS 

(indx,  1) 

= 

hpbwMj  r; 

1732 

hpbwMjrS 

(indx,  2) 

= 

0; 

1733 

hpbwMnrS 

(indx,  1) 

= 

hpbwMnr; 

1734 

hpbwMnrS 

(indx,  2) 

* 

0; 

1735 

rolls 

(indx,  1) 

= 

roll  ; 

1736 

rolls 

(indx,  2) 

= 

0; 

1737 

directivityDBS 

(indx,  1) 

= 

directivityDB; 

1738 

directivityDBS 

(indx,  2) 

= 

0; 

1739 

powe  rMai nVi sbDBS 

(indx,  1) 

= 

powerMainVisbDB 

1740 

powerMainVisbDBS 

(indx,  2) 

= 

0; 

1741 

powerMainSideDBS 

(indx,  1) 

= 

powe rMainSideDB 

1742 

powe  rMainS i deDBS 

(indx,  2) 

= 

0; 

1743 

powerVi sbSideDBS 

(indx,  1) 

= 

powerVisbSideDB 

1744 

powe  rVi sbSideDBS 

(indx,  2} 

= 

0; 

1745 

powerSideAvgDBS 

(indx,  1) 

- 

powers ideAvgDB ; 

1746 

powerSideAvgDBS 

(indx,  2) 

0; 

1747 

slNrstDistS 

(indx,  1) 

** 

slNrstDist; 

1748 

slNrstDistS 

(indx,  2) 

- 

0; 

1749 

slNrstPowrDBS 

(indx,  1) 

“ 

slNrst PowrDB; 

1750 

slNrstPowrDBS 

(indx,  2) 

* 

0; 

1751 

slLgstDistS 

(indx,  1) 

* 

slLgstDist; 

1752 

slLgstDistS 

(indx,  2) 

0; 

1753 

slLgstPowrDBS 

(indx,  1) 

" 

slLgstPowrDB; 

1754 

slLgstPowrDBS 

(indx,  2) 

* 

0; 

1755 

1756 

else 

numAcc  (indx)  =  i 

numAcc  (indx) 

+ 

1; 

1757 

factV  =  (numAcc 

(indx)  -  2)  / 

(numAcc  (indx)  - 

1788 

1789 

1790 

1791 

1792 

1793 

1794 

1795 

1796 

1797 

1798 

1799 

1800 
1801 
1802 

1803 

1804 

1805 

1806 

1807 

1808 

1809 

1810 
1811 
1812 

1813 

1814 

1815 

1816 

1817 

1818 

1819 

1820 


-  1); 


1758 

dev  =  px  -  pxS 

(indx,  1); 

1759 

pxS 

(indx,  2)  =  pxS 

(indx,  2) 

* 

factV  + 

1760 

pxS 

(indx,  1)  =  pxS 

(indx,  1) 

+ 

1761 

dev  =  py  -  pyS 

(indx,  1); 

1762 

pyS 

(indx,  2)  =  pyS 

(indx,  2) 

+ 

factV  + 

1763 

pys 

(indx,  1)  -  pyS 

(indx,  1) 

+ 

1764 

dev  -  pz  -  pzS 

(indx,  1); 

1765 

pzS 

(indx,  2)  *  pzS 

(indx,  2) 

* 

factV  + 

1766 

pzS 

(indx,  1)  *»  pzS 

(indx,  1) 

+ 

1767 

dev  -  errPoint 

-  errPointS  (indx,  1); 

1768 

errPointS 

(indx,  2)  *  errPointS 

(indx,  2) 

* 

factV  + 

1769 

errPointS 

(indx,  1)  *  errPointS 

(indx,  1) 

+ 

1770 

S  dev  -  errPointUnc  -  errPointUncS  (indx,  1); 

1771 

V,  errPointUncS 

(indx,  2)  *  errPointUncS 

(indx,  2) 

factV  + 

1772 

V.  errPointUncS 

(indx,  1)  -  errPointUncS 

(indx,  1) 

+ 

1773 

dev  *  beamPowerDB  -  beamPowerDBS  (indx,  1); 

1774 

beamPowerDBS 

(indx,  2)  -  beamPowerDBS 

(indx,  2) 

* 

factV  + 

1775 

beamPowerDBS 

(indx,  1)  -  beamPowerDBS 

(indx,  1) 

+ 

1776 

dev  =  beamDepthDB  -  beamDepthDBS  (indx,  1) ; 

1777 

beamDepthDBS 

(indx,  2)  “  beamDepthDBS 

(indx,  2) 

* 

factV  + 

1778 

beamDepthDBS 

(indx,  1)  -  beamDepthDBS 

(indx,  1) 

+ 

1779 

dev  =  hpbwMj r 

-  hpbwMjrS  (indx,  1); 

1780 

hpbwM j rS 

(indx,  2)  -  hpbwMjrS 

(indx,  2) 

* 

factV  4 

1781 

hpbwMj rS 

(indx,  1)  -  hpbwMjrS 

(indx,  1) 

4 

1782 

dev  -  hpbwMnr 

-  hpbwMnrS  (indx,  1); 

1783 

hpbwMnrS 

(indx,  2)  =  hpbwMnrS 

(indx,  2) 

* 

factV  4 

1784 

hpbwMnrS 

(indx,  1)  =  hpbwMnrS 

(indx,  1) 

4 

1785 

dev  =  roll  -  rolls  (indx,  1); 

1786 

rolls 

(indx,  2)  =  rolls 

(indx,  2) 

* 

factV  4 

1787 

rolls 

(indx,  1)  *=  rolls 

(indx,  1) 

4 

devA2 

dev 


devA2 

dev 


devA2 

dev 


devA2 

dev 


devA2 

dev 


devA2 

dev 


devA2 

dev 


devA2 

dev 


/  numAcc 
/  numAcc 


/  numAcc 
/  numAcc 


/  numAcc 
/  numAcc 


/  numAcc 
/  numAcc 


/  numAcc 
/  numAcc 


/  numAcc 
/  numAcc 


/  numAcc 
/  numAcc 


/  numAcc 
/  numAcc 


dev  =  directivityDB  -  directivityDBS  (indx,  1); 
directivityDBS  (indx,  2)  -  directivityDBS  (indx, 

directivityDBS  (indx,  1)  -  directivityDBS  (indx, 

dev  -  powerMainVisbDB  -  powerMainVisbDBS  (indx,  1); 
powerMainVisbDBS  (indx,  2)  *  powerMainVisbDBS  (indx, 
powerMainVisbDBS  (indx,  1)  -  powerMainVisbDBS  (indx, 
dev  -  powerMainSideDB  -  powerMainSideDBS  (indx,  1); 
powerMainSideDBS  (indx,  2)  =  powerMainSideDBS  (indx, 
powerMainSideDBS  (indx,  1)  -  powerMainSideDBS  (indx, 
dev  =  powerVisbSideDB  **  powerVisbSideDBS  (indx,  1); 
powerVisbSideDBS  (indx,  2)  «  powerVisbSideDBS  (indx, 
powerVisbSideDBS  (indx,  1)  =  powerVisbSideDBS  (indx, 
dev  «=  powerSideAvgDB  -  powerSideAvgDBS  (indx,  1); 
powers ideAvgDBS  (indx,  2)  -  powerSideAvgDBS  (indx, 

powerSideAvgDBS  (indx,  1)  =  powerSideAvgDBS  (indx, 

slNrstDist  -  slNrstDistS  (indx,  1); 

( indx,  2)  =  slNrstDistS  (indx, 

(indx,  1)  =  slNrstDistS  (indx, 

dev  -  slNrst PowrDB  -  slNrstPowrDBS  (indx,  1); 
slNrstPowrDBS  (indx,  2)  =  slNrstPowrDBS  (indx, 

(indx,  1)  «  slNrstPowrDBS  (indx, 

slLgstDistS  (indx,  1) ; 

(indx,  2)  =  slLgstDistS  (indx, 

(indx,  1)  =  slLgstDistS  (indx, 

dev  =  slLgstPowrDB  -  slLgs tPowrDBS  (indx,  1)  ; 
slLgstPowrDBS  (indx,  2)  =  slLgs tPowrDBS  (indx, 

slLgs tPowrDBS  (indx,  1)  -  slLgstPowrDBS  (indx, 

end 


*  factV  + 


devA2 

dev 


devA2 

dev 


devA2 

dev 


dev 
slNrstDistS 
slNrstDistS 


slNrstPowrDBS 
dev  *  slLgstDist 
slLgstDistS 
slLgstDistS 


devA2 

dev 


devA2 

dev 


*  factV  + 


devA2 

dev 


devA2 

dev 


/  numAcc 
/  numAcc 


/  numAcc 
/  numAcc 


/  numAcc 
/  numAcc 


/  numAcc 
/  numAcc 


/  numAcc 
/  numAcc 


/  numAcc 
/  numAcc 


/  numAcc 
/  numAcc 


devA2 

dev 


devA2 

dev 


devA2 

dev 


devA2 

dev 


/  numAcc 
/  numAcc 


/  numAcc 
/  numAcc 


/  numAcc 
/  numAcc 


/  numAcc 
/  numAcc 


eval  ((‘save  case*  sprintf  ( *  V.O.Of  * ,  indx) 
end 

clear  factV  dev 


sprintf  ( *9.0  .Of ' ,  rlzNum)  ’  excgSqr*]) 


(indx) ; 
(indx)  ; 

(indx)  ; 
(indx)  ; 

(indx)  ; 
(indx)  ; 

(indx)  ; 
(indx)  ; 

(indx) ; 
(indx)  ; 

(indx) ; 
(indx) ; 

(indx) ; 
(indx) ; 

(indx) ; 
(indx) ; 

(indx) ; 
(indx) ; 

(indx) ; 
(indx) ; 

(indx) ; 
(indx) ; 

(indx); 
(indx) ; 

(indx) ; 
(indx) ; 

(indx) ; 
(indx)  ; 

(indx) ; 
(indx) ; 

(indx) ; 
(indx) ; 

(indx) ; 
(indx) ; 

(indx) ; 
(indx) ; 

(indx) ; 
(indx) ; 
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7  Calculate  spherical  coordinates  of  average  pointing  vector 

7.  Let  <px>,  <py>,  and  <pz>  denote  the  average  direction  cosines  of 
V.  the  pointing  vectors,  and  let  vx,  vy,  and  vz  denote  the 
7,  corresponding  variances.  We  wish  to  express  the  direction  of  the 
V.  average  pointing  vector  |<px>  <py>  <pz>)  in  spherical  coordinates. 

|  First,  note  that  the  average  pointing  vector  has  the  norm 

V.  2  2  2 

7.  p  =  ( <px >  +  <py>  +  <pz>  )  , 

7,  which  is  less  than  one  if  not  all  realizations  are  colinear.  The 
7  spherical  angles  are  then  given  by 

7.  12  2  1/2 

7,  sin  pointPolar  *  -  (<px>  +  <py>  )  and 

7  P 

v,  <py> 

V,  tan  pointAzmth  «  — --  , 

7.  <px> 

V. 

7.  where  we  use  sin  pointPolar  instead  of  the  cosine  for  accuracy 
7.  near  boresight. 

% 

7  Also,  we  wish  to  estimate  the  rms  angular  deviation  of  a 
7  realization  of  the  pointing  vector  from  the  mean.  Adopting  an 
7  unsophisticated  method,  we  add  the  variances  vx  and  vy  to  obtain 
7  an  equivalent  area  in  the  x-y  direction  cosine  plane,  then  divide 
7  the  area  by  cos  PointPolar  to  yield  a  solid  angle  on  the 
7  hemisphere.  The  square  root  of  that  area  yields  the  rms  angular 
7  deviation. 

7 

radSqr  -  pxS  (indx,  1)A2  +  pyS  (indx,  1)A2; 

pointPolar  -  asin  (sqrt  (radSqr  /  (radSqr  +  pzS  {indx,  1)A2))); 
pointAzmth  -  atan2  (pyS  (indx,  1),  pxS  (indx,  1)); 
if  isnan  (pointPolar) 
pointStdDev  -  nan; 
else 

pointStdDev  -  sqrt  (<pxS  (indx,  2)  +  pyS  (indx,  2))  /  cos  (pointPolar)); 
end 

clear  radSqr 

7  Print  performance  characteristics 
7 

7  Generally,  the  following  statements  print  the  results  of  the 
7  analysis  followed  by  the  standard  deviations  of  each  result  in 
7  curly  brackets. 

7 

if  1  7  INPUT  0  to  suppress  output,  1  to  print 

fprintf  (1,  ' \nMeans  and  [std  devs)  for  70. Of  of  70. Of  realizations\n* ,  rlzNum,  numRlz) ; 
fprintf  (1,  ‘beam  direction  ;  (70. 3f,  70. 3 f )  deg,  std  dev  70. 3f  deg\n',  ... 

pointPolar  /  rpd,  (pointAzmth  +  azmthOf fst)  /  rpd,  pointStdDev  /  rpd) ; 
fprintf  (1,  ’pointing  error  :  78. 4 f  (70.4fJ  deg\n’,  — 

errPointS  (indx,  1)  /  rpd,  sqrt  (errPointS  (indx,  2))  /  rpd); 

7  fprintf  (1,  'pntng  error  unc:  78. 4f  (70.4fJ  deg  (2  sigma) \n*,  ... 

7  2  *  errPointUncS  (indx,  1)  /  rpd,  2  *  sqrt  (errPointUncS  (indx,  2))  /  rpd); 

fprintf  (1,  'peak  power  dens:  7.7. 3  f  (70.3f  )  dB\n*,  ... 

beamPowerDBS  (indx,  1),  sqrt  (beamPowerDBS  (indx,  2))); 
fprintf  (1,  ’beam  depth  :  76. 2f  (70. 2f  )  dB  re  peak\n\  ... 

beamDepthDBS  (indx,  1),  sqrt  (beamDepthDBS  (indx,  2))); 
fprintf  (1,  ’beam  width  :  (76. 3f  (70. 3f  ),  70. 3f  [70. 3f ] )  deg\n’,  ... 

hpbwMjrS  (indx,  1)  /  rpd,  sqrt  (hpbwMjrS  (indx,  2))  /  rpd,  ... 
hpbwMnrS  (indx,  1)  /  rpd,  sqrt  (hpbwMnrS  (indx,  2))  /  rpd); 
fprintf  (1,  'beam  roll  :  76. 2f  (70. 2f  )  deg\n*,  ... 

(rolls  (indx,  1)  +  azmthOffst)  /  rpd,  sqrt  (rolls  (indx,  2))  /  rpd); 
fprintf  (1,  ’directivity  :  77. 3f  (70. 3f  ]  dB\n',  ... 

direct ivityDBS  (indx,  1),  sqrt  (directivityDBS  (indx,  2))); 
fprintf  (1,  'power  ratio  m/v:  77. 3f  (70. 3f  )  dB\n’,  ... 

powerMainVisbDBS  (indx,  1),  sqrt  (powerMainVisbDBS  (indx,  2))); 
fprintf  (1,  'power  ratio  m/s:  77. 3f  (70. 3f  ]  dB\n*,  ... 

powerMa inSide DBS  (indx,  1),  sqrt  (powerMainSideDBS  (indx,  2))); 
fprintf  (1,  'power  ratio  v/s:  77. 3f  (70. 3f  ]  dB\n’,  ... 

powerVisbSideDBS  (indx,  1),  sqrt  (powerVisbSideDBS  (indx,  2))); 
fprintf  (1,  'avg  sidelobe  :  76. 2f  (70. 2f  ]  dB  re  peak\n',  ... 

powers ideAvgDBS  (indx,  1),  sqrt  (powerSideAvgDBS  (indx,  2))); 
f print f  (1,  ' nrst  sidelobe  :  76. 2f  (70. 2f]  dB  re  peak,  70. 2f  (70. 2f)  deg  off  beam\n',  ... 

SlNrstPowrDBS  (indx,  1),  sqrt  (slNrstPowrDBS  (indx,  2}),  ... 

slNrstDistS  (indx,  1)  /  rpd,  sqrt  (slNrstDistS  (indx,  2))  /  rpd); 
fprintf  (1,  ’ lgst  sidelobe  :  76. 2f  (70. 2f)  dB  re  peak,  70. 2f  (70. 2f]  deg  off  beam\n',  ... 

s 1 Lgs t PowrDBS  (indx,  1),  sqrt  (slLgstPowrDBS  (indx,  2)),  ... 
slLgstDistS  (indx,  1)  /  rpd,  sqrt  (slLgstDistS  (indx,  2))  /  rpd); 

end 

end  7  loop  over  realizations 

7  Plot  performance  characteristics  as  function  of  independent  variable 
7  (summary  plot) . 

7 

7  If  more  than  one  realization  has  been  accumulated  for  any  value  of 
7  the  independent  variable,  the  mean  is  plotted  with  uncertainty  bars. 

7  The  extension  of  the  uncertainty  bar  above  the  mean  equals  one 
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V.  standard  deviation,  and  likewise  below  the  mean.  Otherwise,  only 
V,  the  values  for  the  one  realization  are  plotted. 


V.  Plots  that  show  more  than  one  measure  distinguish  them  by  color  or 
*„  line  style  and  may  use  both  a  left  and  right  axis.  The  color  or 
style  and  axis  for  each  measure  is  given  in  codes  in  parentheses  in 
V.  the  title  of  the  plot.  The  first  code  abbreviates  the  color  or  line 
style: 


R  red  -  solid 

V.  G  green  --  dashed 

■A  B  blue  :  dotted 

V.  C  cyan 

V.  M  magenta 

Y  yellow 
V.  K  black 

V,  W  white 


l  and  the  second  letter  indicates  the  axis,  L  for  left  and  R  for 
V.  right. 


if  indVarLen  >  1 
figure  (figSum); 
elf ; 

axesSum  -  zeros  (12,  1);  Vi  space  for  axes  handles 

I  Axes  1:  pointing  error 

% 

subplot  {4,  2,  1); 
axesSum  (1)  =  gca; 

if  any  (numAcc  >1)  , 

hline  -  errorbar  (indVar,  errPointS  {:,  1)  /  rpd,  sqrt  (errPointS  2))  /  rpd); 

set  (hline,  discrim,  discrimValue  {1});  %  possibly  override  errorbar ’s  default  solid  line  style 

set  (hline  (1),  ’linestyle’,  •-•);  *  but  leave  the  error  bars  themselves  solid 

else 

plot  (indVar,  errPointS  (:,  1)  /  rpd) ; 
end 

ylabel  ( * (deg) * ) ; 
title  {‘pointing  error’) ? 

‘A  Axes  2  and  3:  beam  power  and  directivity 


subplot  (4,  2,  2); 
axesSum  (2)  *=  gca; 
if  any  (numAcc  >1) 

hline  -  errorbar  (indVar,  beamPowerDBS  (:,  1),  sqrt  (beamPowerDBS  (:,  2))); 
set  (hline, '  discrim,  discrimValue  ID); 
set  (hline  (1),  ’linestyle*, 
else 

plot  (indVar,  beamPowerDBS  (:,  1)); 
end 

ylabel  (*(dB  re  coherent)*); 

title  (strcat  {’peak  power  dens  {’,  discrimName  |1),  ... 

•  l) ,  directivity  (’,  discrimName  (2),  ’  R)*)); 
axesSum  (3)  *  axes  (’position*,  get  (gca,  ’position’)); 

if  any  (numAcc  >1)  .  .  _ _  ,  .... 

hline  -  errorbar  (indVar,  directivityDBS  (:,  1),  sqrt  (di recti vityDBS  {:,  2))); 
set  (hline,  discrim,  discrimValue  (2)); 
set  (hline  (1),  ’linestyle’,  *-’); 
else 

plot  (indVar,  directivityDBS  (:,  1),  discrim,  discrimValue  (2>); 
end 

set  (gca,  ’yAxisLocation’ ,  ’right*,  ’color*,  ’none’); 


% 

'i  Axes  4  and  5:  beam  widths 
% 

subplot  (4,  2,  3); 
axesSum  (4)  =  gca; 
if  any  (numAcc  >  1) 

hline  =  errorbar  (indVar,  hpbwMjrS  (:,  1)  /  rpd,  sqrt  (hpbwMjrS  {:, 

set  (hline,  discrim,  discrimValue  (D); 
set  (hline  (1),  'linestyle’,  ’-*); 
else 

plot  (indVar,  hpbwMjrS  (:,  1)  /  rpd); 
end 

ylabel  ( * (deg) ’ ) ; 

title  (strcat  (’hpbw  major  (’,  discrimName  (1),  — 

«  L)  &  minor  {’,  discrimName  (2),  ’  R) ’ ) ) ; 

axesSum  (5)  -  axes  {’position*,  get  (gca,  ’position*)); 
if  any  (numAcc  >  1) 

hline  =  errorbar  (indVar,  hpbwMnrS  (:,  1)  /  rpd,  sqrt  (hpbwMnrS  (:, 

set  (hline,  discrim,  discrimValue  (2)); 
set  (hline  (1),  ’linestyle’,  ’-’); 
else 

plot  (indVar,  hpbwMnrS  (:,  1)  /  rpd,  discrim,  discrimValue  (2)); 


2)}  /  rpd}; 


2))  /  rpd) ; 


end 

set  (gca,  ’yAxisLocation*,  ’right’,  ’color’,  ’none’); 
yLimL  =  get  (axesSum  (4),  ’yliro’); 
yLimR  «=  get  (axesSum  (5),  ’ylim’); 
if  yLimL  (1)  <  yLimR  (2) 

set  (axesSum  (4),  ’xlimmode’,  ’manual*,  ’ylim’,  [yLimR(l)  yLimL(2)]); 
set  (axesSum  (5),  ’xlimmode’,  ’manual’,  ’ylim’,  (yLimR(l)  yLimL(2)]); 


,9  c. 8' 

Ty'.P, 
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V.  Setting  xlimraode  to  manual  prevents  rescaling  of  the  x  axis  when 
V.  the  y  axis  is  changed, 
end 

",  Axes  6:  beam  roll 

subplot  ( 4 ,  2,  4); 
axesSum  (6)  =  gca; 
if  any  (numAcc  >  1) 

hline  =  errorbar  (indVar,  rolls  (:,  1}  /  rpd,  sqrt  (rolls  (:,  2))  /  rpd)  ; 

set  (hline,  discrim,  discrimValue  (1}); 
set  (hline  (1),  'linestyle',  '-*); 
else 

plot  (indVar,  rolls  (:,  1)  /  rpd); 
end 

ylabel  ( ’ (deg) * ) ; 
title  ('beam  roll'); 

ft 

V.  Axes  7:  beam  depth 
ft 

subplot  (4,  2,  5); 
axesSum  (7)  «  gca; 
if  any  (numAcc  >  1) 

hline  -  errorbar  (indVar,  beamDepthDBS  (:,  1),  sqrt  (beamDepthDBS  (:,  2))); 
set  (hline,  discrim,  discrimValue  (1>); 
set  (hline  (1),  'linestyle*,  *-'); 
else 

plot  (indVar,  beamDepthDBS  (:,  1)); 
end 

ylabel  ( * (dB  re  peak)*); 
title  (’beam  depth’); 
ft 

ft  Axes  8  and  9:  power  ratios 
ft 

subplot  (4,  2,  6); 
axesSum  (8)  «  gca; 
if  any  (numAcc  >  1) 

hline  -  errorbar  (indVar,  powerMainVisbDBS ( : , 1) ,  sqrt  (powerMainVisbDBS (:, 2) )) ; 
set  (hline,  discrim,  discrimValue  {1}); 
set  (hline  (1),  'linestyle*,  '-*); 
else 

plot  (indVar,  powerMainVisbDBS  (:,  1) )  ; 
end 

ylabel  C(dB)'); 

title  (strcat  (’power  m/v  (*,  discrimName  (1),  ... 

*  L),  m/s  (',  discrimName  (2),  ... 

•  R),  v/s  (',  discrimName  (3),  '  R) ' ) ) ; 

axesSum  (9)  -  axes  ('position',  get  (gca,  'position')); 

if  any  (numAcc  >  1) 

hline  -  errorbar  (indVar,  powerMainSideDBS ( : ,  1) ,  - 

sqrt  (powerMainSideDBS (:  ,2) )) ; 
set  (hline,  discrim,  discrimValue  {2}); 
set  (hline  (1),  'linestyle', 
set  (gca,  'nextplot',  'add'); 

hline  -  errorbar  (indVar,  powerVisbSideDBS ( : ,  1) ,  ... 

sqrt  (powerVisbSideDBS (:  ,2) ))  ; 
set  (hline,  discrim,  discrimValue  (3)); 
set  (hline  (1),  'linestyle',  '-'); 
set  (gca,  'nextplot*,  'replace'); 
else 

plot  (indVar,  powerMainSideDBS (;, 1) ,  discrim,  discrimValue  (2>); 
set  (gca,  'nextplot*,  'add'); 

plot  (indVar,  powerVisbSideDBS (:, 1) ,  discrim,  discrimValue  (3)); 
set  (gca,  'nextplot',  'replace'); 
end 

set  (gca,  'yAxisLocation* ,  'right*,  'color',  'none'); 
ft 

ft  Axes  10  and  11:  sidelobe  powers 
ft 

subplot  (4,  2,  7); 
axesSum  (10)  -  gca; 
if  any  (numAcc  >  1) 

hline  -  errorbar  (indVar,  slLgstPowrDBS  (:,  1),  sqrt  (slLgstPowrDBS  (:,  2))); 
set  (hline,  discrim,  discrimValue  (1)); 
set  (hline  (1),  'linestyle', 
set  (gca,  'nextplot',  'add'); 

hline  =  errorbar  (indVar,  slNrstPowrDBS  (:,  1),  sqrt  (slNrs tPowrDBS  (:,  2))); 
set  (hline,  discrim,  discrimValue  (2)); 
set  (hline  (1),  'linestyle',  '-'); 
set  (gca,  'nextplot',  'replace'); 
else 

plot  (indVar  *  ones  (1,  2),  ... 

(slLgstPowrDBS ( : , 1)  slNrstPowrDBS ( : , 1) ) ) ; 

end 

ylabel  (' (dB  re  peak) ’) ; 

title  (strcat  ('sidelobe  power:  lgst  (',  discrimName  (1),  ... 

*  L) ,  nrst  (’,  discrimName  (2),  ... 

*  L) ,  avg  (’,  discrimName  (3),  *  R)’)); 

axesSum  (11)  -  axes  ('position',  get  (gca,  'position')); 

if  any  (numAcc  >  1) 

hline  «=  errorbar  (indVar,  powerSideAvgDBS  (:,  1),  sqrt  (powerSideAvgDBS  (:,  2))); 
set  (hline,  discrim,  discrimValue  (3)); 
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set  (hline  (1),  'linestyle', 
else 

plot  (indVar,  powers ideAvgDBS  1),  discrim,  discrimValue  (3)); 

end 

set  (gca,  1 yAxisLocation * ,  ‘right*,  'color',  'none'); 


ft  Axes  12:  sidelobe  distances 
ft 

subplot  (4,  2,  8) ; 
axesSum  (12)  =  gca; 
if  any  (numAcc  >  1) 

hline  =  errorbar  (indVar,  slLgstDistS  (:,  1)  /  rpd,  sqrt  (slLgstDistS  {:, 
set  (hline,  discrim,  discrimValue  {1}); 
set  (hline  (1),  'linestyle',  '-'); 
set  (gca,  'nextplot',  'add'); 

hiine  -  errorbar  (indVar,  slNrstDistS  (:,  1)  /  rpd,  sqrt  (slNrstDistS  (:, 
set  (hline,  discrim,  discrimValue  {2}); 
set  (hline  (1),  'linestyle',  '-'); 
set  (gca,  'nextplot',  'replace'); 
else 

plot  (indVar  *  ones  (1,  2),  ... 

[slLgstDistS (:,1)  slNrstDistS (:, 1) ]  /  rpd); 

end 

ylabel  (*{deg)'); 

title  (strcat  ('sidelobe  distance:  Igst  (*,  discnmName  (1),  ... 

'),  nrst  (*,  discrimName  (2),  *)')); 


2))  /  rpd) 


2) )  /  rpd) 


% 

ft  Touch  up 
ft 

set  (f indob j  (gcf,  'type',  'axes'),  'xlim',  ... 
(min (indVar)  max (indVar) ]  ... 

+  0.1  *  (max  (indVar)  -  min  (indVar))  *  [-1  U); 
axesSumTitle  *=  axes  ('position',  [0011],  ... 

'color',  'none',  'visible',  'off*,  ... 

' def aultTextFontSize ' ,  10,  ... 

•def aultTextHorizontalAlignment* ,  'center* ) ; 
text  (0.5,  0.0S,  indVarName,  ... 
'horizontalAlignment 'center'); 

end 

clear  hline 


ft  make  x-axis  limits  uniform 


ft  create  invisible  axes 
ft  for  titling 

ft  display  name  of  independent 
ft  variable 


ft  Plot  last  realization 
ft 

ft  Three  projections  of  the  hemisphere  are  available:  Lambert, 
ft  stereographic,  and  orthographic, 
ft 

ft  Lambert  projection: 

ft  . 

ft  The  Lambert  projection  preserves  the  relative  areas  of  portions  of 
ft  the  hemisphere.  That  is,  the  ratio  of  areas  of  two  regions  on  the 
ft  projection  is  the  same  as  on  the  hemisphere.  The  azimuth  coordinate 
ft  of  a  point  in  the  projection  is  the  same  as  its  azimuth  coordinate 
ft  on  the  hemisphere,  while  its  radius  r  from  the  center  of  the 
ft  projection  is  related  to  the  polar  angle.  This  relationship  may  be 
ft  derived  by  setting  the  spherical  surface  area  sin  polar  d (polar) 
ft  d(azmth)  equal  to  a  constant  times  the  planar  surface  area  r  dr 
ft  d(azmth)  and  integrating.  The  radius  is  then  given  by 


ft  1/2  polar 

ft  r  =  2  sin - . 

ft  2 

ft  For  computer  graphics  the  Cartesian  coordinates  are  more  convenient; 
ft  they  are 
ft 

ft  u  »  r  cos  azmth  *  R  cx 

ft 

ft  v  «  r  sin  azmth  =  R  cy 


ft  where 

ft 

ft  -1/2  polar  -1/2 

ft  R  =  2  sec - -  (1  +  cz) 

ft  2 

ft  For  points  outside  visible  space  (that  is,  for  cx^2  +  cy  2  >  1),  cz 
ft  is  zero  so  that  R  =  1  and  no  transformation  is  applied. 

ft  Stereographic  projection: 

V.  The  stereographic  projection  preserves  angles  on  the  hemisphere.  To 
ft  derive  the  governing  relationship  for  the  projection,  use  the  same 
ft  azimuthal  angle  for  the  projection  as  for  the  hemisphere,  and  let 
ft  the  radius  be  a  function  of  the  polar  angle:  r  =  f  (polar).  Equate 
ft  the  aspect  ratios  of  orthogonal  derivatives,  as 
ft 

V.  d  (polar)  dr 

ft  - “ - 

ft  sin  polar  d (azmth)  r  d (azmth) 
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V.  f  *  (polar)  d(polar) 


V  f  (polar)  d(azmth) 

(where  f*  is  the  first  derivative),  rearrange,  and  integrate  to 
v.  obtain 

polar  1  -  cos  polar 

f  (polar)  =  tan - * - , 

2  sin  polar 

¥.  up  to  an  arbitrary  multiplicative  constant.  The  projected  Cartesian 
V.  coordinates  of  a  point  (cx,  cy,  cz)  on  the  sphere  are 

V  u  *  r  cos  azmth  =  R  cx 


V.  v  *  r  sin  azmth  *=  R  cy 

9. 

V,  where 


V.  f  (polar)  1  1 

V,  R  - - - - - - - 

V,  sin  polar  1  +  cos  polar  1  +  cz 

V, 

%  For  points  outside  visible  space  (that  is,  for  cxA2  +  cyA2  >1),  cz 
%  is  zero  so  that  R  -  1  and  no  transformation  is  applied. 

9. 

9,  The  center  of  projection  is  opposite  boresight  (cx  *  cy  -  0,  cz  « 

V,  -1),  and  with  the  above  choice  of  multiplicative  constant,  the  plane 
V.  of  projection  is  the  cz  -  0  plane. 


V.  Orthographic  projection: 

9. 

9.  The  orthographic  projection  gives  a  3D  view  of  the  hemisphere. 
4 

if  1  %  INPUT 

INPUT 


flat* 

0; 


proj  -  1; 


coordSys  ■  1; 

faceColor  - 
pointZoom  * 

9.  The  "show** 
showBeamRegion  “ 
showWidthRegion  * 
showWidthContAct  - 
showWidthContFit  • 
showPointGrid  - 
showPoint  « 

showPointUnc  » 
showSidelobeGrid  - 
showSidelobeNrst  - 
showSidelobeLgst  * 


%  INPUT  1  to  plot,  0  to  skip 
4  INPUT  1  for  2D  equal-area  Lambert 

9.  2  for  2D  stereographic 

4  3  for  orthographic  (3D  hemisphere) 

9.  INPUT  1  for  grid  of  spherical  coordinates 

9  2  for  grid  of  traditional  coordinates 

9  INPUT  'flat'  or  'interp*  shading 

9  INPUT  magnification  factor  or  0  for  centered  full  view 


inputs  below  are  coded  only  for  2D  views. 


9  INPUT  1  to  show  main  beam  region 


region  above  width  contour 

actual  width  contour 

fitted  width  contour  (ellipse) 

main  beam  peak  on  sampled  grid 

interpolated  main  beam  peak  (pointing  vector) 

axes  of  uncertainty  ellipse  of  pointing  vector 

sidelobe  peaks  on  sampled  grid 

nearest  sidelobe 

largest  sidelobe 


9,  Prepare  for  plotting 

9. 

if  strcmp  (faceColor,  'interp’) 

di  rCosXMtxSurf  *=  dirCosXMtx;  9.  use  true  grid 

dirCosYMtxSurf  “  dirCosYMtx; 
di rCosZMtxSurf  -  dirCosZMtx; 
visBoolSurf  -  visBool  ; 
elseif  strcmp  (faceColor,  ’flat’) 

di rCosXMtxSurf  -  ones  (ay,  1)  *  dirCosShiftX;  *  use  shifted  grid  to  center 
dirCosYMtxSurf  -  dirCosShiftY  *  ones  (1,  ax);  9.  patches  on  data  points 
radSqr  -  dirCosXMtxSurf . A2  +■  dirCosYMtxSurf . A2; 
visBoolSurf  -  (radSqr  <  1); 
di rCosZMtxSurf  -  zeros  (ay,  ax); 

di rCosZMtxSurf  (visBoolSurf)  -  sqrt  (1  -  radSqr  (visBoolSurf)); 
clear  radSqr; 
else 

error  ('Illegal  value  of  faceColor.'); 


end 

if  strcmp  (faceColor,  'interp* 
gBlnklndx  «  zeros  (ay,  ax); 
gBlnklndx  (l:ay-l,  l:ax-l)  ■ 
gBlnklndx  (l:ay-l,  2:ax  )  = 
gBlnklndx  (2:ay  ,  l:ax-l)  = 

gBlnklndx  (2:ay  ,  2 : ax  )  - 

else 

gBlnklndx  *  visBoolSurf; 
gBlnklndx  (l:ay-l,  l:ax  )  * 


(l:ay-l,  l:ax-l) 
1 :  ax-1 ) 


gBlnklndx 
gBlnklndx  (l:ay 
end 

gBlnklndx  =  -gBlnklndx; 
gZeroIndx  •*  (gSqr  ««  0); 
gOKIndx  -  find  (-(gBlnklndx  | 
gZeroIndx  -  find  (gZeroIndx); 
gBlnklndx  -  find  (gBlnklndx); 
if  1 


)  ",  identify  visible  points  plus  those  immediately 

9.  or  diagonally  adjacent  (Boolean  for  now) 
visBool  (2:ay,  2:ax); 

gBlnklndx  (l:ay-l,  2:ax  )  |  visBool  (2:ay  ,  l:ax-l); 

(2:ay  ,  l:ax-l)  J  visBool  (l:ay-l,  2:ax  ); 

(2:ay  ,  2:ax  )  l  visBool  (l:ay-l,  l:ax-l); 

9.  identify  visible  points  plus  those  immediately 
9.  but  not  diagonally  adjacent  (Boolean  for  now) 
(l:ay-l,  l:ax  )  I  visBoolSurf  (2:ay  ,  l:ax); 
<l:ay-l,  l:ax-l)  1  visBoolSurf  (2:ay  ,  2:ax); 

(l:ay  ,  l:ax-l)  I  visBoolSurf  (l:ay  ,  2:ax); 


gBlnklndx 

gBlnklndx 


gBlnklndx 

gBlnklndx 

gBlnklndx 


9.  0  in  visible  region  plus  appropriate  border 
9,  1  where  zero,  0  elsewhere  (Boolean  for  now) 
gZeroIndx));  9,  1  for  nontrivial  values 
9.  convert  to  indices 


V.  Decibel  plot 


*  INPUT  0  for  linear,  1  for  decibel  plot 
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V.  INPUT  dB  range 

x  {dim  (1),  10  *  log  10  (gSqr  (gOKIndx))); 

dim  (1)  *  ones  (size  (gZeroIndx));  %  condition  log  of  0 


Set  appropriate  value  for  invisible  points 
ones  (size  (gBlnklndx)); 


gPlot  -  zeros  (ay,  ax); 
dim  =  (-50  0]; 

gPlot  (gOKIndx)  =  max  (dim  (1),  10  *  loglO 
gPlot  (gZeroIndx)  =  dim  (1)  *  ones  (size  (< 
else 

v  Linear  plot 
gPlot  =  gSqr; 
dim  =  (0  1] ; 
end 

if  invertBkgd 

gPlot  (gBlnklndx)  =  dim  (2)  *  ones  (size  (< 
else 

gPlot  (gBlnklndx)  =  dim  (1)  *  ones  (size  (gBlnklndx)); 
end 

dear  visBoolSurf  gBlnklndx  gZeroIndx  gOKIndx 
figure  (figPat) ; 

fore  -  get  (figPat,  ’defaultLineColor ' ) ;  * 

if  invertBkgd  ‘i  choose  fore-  and  background  < 
fore  -  1  -  fore; 
end 

back  -  1  -  fore; 
switch  coordSys 
case  1 

gridlPolar  -  (0  ;  10  :  90)  *  rpd;  V, 

gridlAzmth  -  (0  :  5  :  360)  *  rpd  -  azmth 

grid2Polar  -  (0  :  5  :  90)  *  rpd;  % 

grid2Azmth  -  (0  :  30  :  360)  *  rpd  -  azmth 


neColor*);  V,  expect  black  or  white 

and  background  colors  compatible  with  color  map 


10  :  90) 

5  :  360) 
5  :  90) 

30  :  360) 


rpd;  %  lines  of  constant  polar  angle 

rpd  -  azmthOffst; 

rpd;  %  lines  of  constant  azimuth 

rpd  -  azmthOffst; 


(gridlPolar,  gridlAzmth]  *  meshgrid  (gridlPolar,  gridlAzmth) ; 
farid2 Azmth,  grid2Polar]  =meshgrid  (grid2Azmth,  grid2Polar) ; 


[grid2 Azmth,  grid2Pol 
case  2 


gridlAz  -  (-90  :  10  :  90)  *  rpd;  *  lines  o 

gridlEl  -  (-90  :  5  ;  90)  *  rpd; 

grid2Az  =  (-90  :  10  ;  90)  *  rpd;  -  lines  o 

grid2El  -  (-90  :  10  ;  90)  *  rpd; 

(gridlAz,  gridlEl)  -  meshgrid  (gridlAz,  gridlEl); 
(grid2El,  grid2Az)  -  meshgrid  (grid2El,  grid2Az); 
gridlPolar  *  acos  (cos  (gridlEl)  .*  cos  (gridlAz 

gridlAzmth  ■*  atan2  (tan  (gridlEl),  -sin  (gridlAz 

grid2Polar  -  acos  (cos  <grid2El)  .*  cos  (grid2Az 

grid2Azmth  -  atan2  (tan  (grid2El),  -sin  (grid2Az 

otherwise 

error  (‘Invalid  value  of  coordSys.*); 


%  lines  of  constant  azimuth 


lines  of  constant  elevation 


*  cos  (gridlAz)); 

-sin  (gridlAz))  -  azmthOffst; 

*  cos  (grid2Az) ) ; 

-sin  (grid2Az) )  -  azmthOffst; 


V.  Plot  array  factor 


if  proj  ==*  1 

% 

l  Lambert  projection 

radFactSurf  -  1  ./  sqrt  (1  +  dirCosZMtxSurf ) ;  %  radius  factor  for  surface  plot 

hSurf  «  surf  (radFactSurf  .*  dirCosXMtxSurf ,  ... 

radFactSurf  .*  dirCosYMtxSurf ,  gPlot); 
gridlRadFact  =  sqrt  (2)  *  sin  (gridlPolar  /  2); 
line  (gridlRadFact  .*  cos  (gridlAzmth),  ... 

gridlRadFact  .4  sin  (gridlAzmth),  ... 

dim  (2)  *  ones  (size  (gridlRadFact)),  ‘color’,  fore); 
grid2RadFact  =  sqrt  (2)  *  sin  (grid2Polar  /  2); 
line  (grid2RadFact  .*  cos  (grid2Azmth) ,  ... 

grid2RadFact  .*  sin  (grid2Azmth) ,  ... 

dim  (2)  *  ones  (size  (grid2RadFact ) ) ,  'color*,  fore); 
set  (gca,  • drawmode-,  'fast');  *  no  hidden  objects  to  worry  about 

xylim  -  (-1  1] ; 
zlim  *=  dim; 

vi pwAz  “0*  %  Matlab  azimuth  and 

viewEl  -  90  *  rpd;  *  elevation  (but  in  radians) 

elseif  proj  —  2 


%  Matlab  azimuth  and 
V.  elevation  (but  in  radians) 


color',  fore); 


Vi  Stereographic  projection 

radFactSurf  =  1  ./  (1  +  dirCosZMtxSurf);  V,  radius  factor  for  surface  plot 

hSurf  =  surf  (radFactSurf  .*  dirCosXMtxSurf,  ... 

radFactSurf  .*  dirCosYMtxSurf,  gPlot); 
gridlRadFact  =  tan  (gridlPolar  /  2); 
line  (gridlRadFact  .*  cos  (gridlAzmth),  ... 

gridlRadFact  .*  sin  (gridlAzmth),  ... 

dim  (2)  *  ones  (size  (gridlRadFact)),  'color',  fore); 
grid2RadFact  -  tan  (grid2Polar  /  2) ; 
line  (grid2RadFact  .*  cos  (grid2Azmth) ,  ... 

grid2RadFact  . *  sin  (grid2Azmth) ,  ... 

dim  (2)  *  ones  (size  (grid2RadFact) ) ,  'color' ,  fore); 
set  (gca,  'drawmode',  'fast');  %  no  hidden  objects  to 

xylim  =  (-1  1] ; 
zlim  =  dim; 

viewAz  -  0;  Matlab  azimuth  and 

viewEl  -  90  *  rpd;  *  elevation  (but  in  r 

elseif  proj  ==  3 

l 

V.  Orthographic  projection  (3D  hemisphere) 

float  “  1.02;  «  radius  of  annotations  relative  to  unit  hemisphere 

hSurf  =  surf  (dirCosXMtxSurf,  dirCosYMtxSurf,  dirCosZMtxSurf,  gPlot); 


hidden  objects  to  worry  about 


V.  Matlab  azimuth  and 
V.  elevation  (but  in  radians) 
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gridlRadFact  -  sin  (gridlPolar) ; 

line  (gridlRadFact  .*  cos  (gridlAzmth}  *  float,  ... 

gridlRadFact  .*  sin  (gridlAzmth)  *  float,  ... 

cos  (gridlPolar)  *  float,  ’color*,  fore); 
grid2RadFact  *  sin  (grid2Polar) ; 

line  (grid2RadFact  .*  cos  (grid2Azmth)  *  float,  ... 

grid2RadFact  .*  sin  (grid2Azmth)  *  float,  ... 

cos  (grid2Polar)  *  float,  'color',  fore); 
set  (gca,  'drawmode',  'normal');  V.  remove  hidden  objects 

xylim  =  (-1  1]  *  float; 
zlim  =  (-1  1)  *  float; 

switch  1  V.  INPUT  1  to  look  down  main  beam;  2,  down  boresight;  3,  custom 
case  1 

'}.  Look  down  main  beam 

viewAz  =  pi/2  +  steerAzmth  +  azmthOffst;  ?.  Matlab  azimuth  and 

viewEl  =  pi/2  -  steerPolar;  V.  elevation  (but  in  radians) 

case  2 

V.  Look  down  boresight 

viewAz  «=  0;  V.  Matlab  azimuth  and 

viewEl  ■  90  *  rpd;  V,  elevation  (but  in  radians) 

set  (gca,  'drawmode',  'fast');  '*  no  hidden  surfaces 

otherwise 

7,  Look  somewhere 

viewAz  -  60  *  rpd;  V.  INPUT  custom  view  azimuth 

viewEl  -  30  *  rpd;  *  and  elevation  (Matlab  coordinates) 

end 

end  7.  plotting 

Annotate  plot  (2D  only) 

if  (proj  —  1)  I  (proj  «  2) 

switch  proj  V.  set  appropriate  radius  factor 
case  1 

radFact  «  1  ./  sqrt  (1  +  dirCosZMtx); 
case  2 

radFact  -  1  ./  (1  +  dirCosZMtx); 

end 

if  showBeamRegion  &  beamExist 

line  (radFact  (beamlndx)  .*  dirCosXMtx  (beamlndx),  ... 
radFact  (beamlndx)  .*  dirCosYMtx  (beamlndx),  ... 
dim  (2)  *  ones  (size  (beamlndx)),  ... 

'linestyle',  'none*,  'marker',  'color',  fore); 

end 

if  showWidthRegion  &  capClosed 

line  (radFact  (caplndx)  .*  dirCosXMtx  (caplndx),  ... 

radFact  (caplndx)  . *  dirCosYMtx  (caplndx),  ... 
dim  (2)  *  ones  (size  (caplndx)),  ... 

'linestyle*,  'none',  'marker*,  'x',  'color',  fore); 

end 

if  showWidthContAct  &  capClosed 
switch  proj 
case  1 

radFactCapAct  -  1  ./  sqrt  (1  +  capContZ); 
case  2 

radFactCapAct  -  1  ./  (1  +  capContZ) ; 

end 

line  (radFactCapAct  .*  capContX,  radFactCapAct  .*  capContY,  ... 

dim  (2)  *  ones  (size  (capContX)),  'linestyle',  ';*,  'color',  back); 
clear  radFactCapAct 
end 

if  showWidthContFit  t  capClosed 
switch  proj 
case  1 

radFactCapFit  -  1  ./  sqrt  (1  +  cFitZ) ; 
case  2 

radFactCapFit  »  1  ./  (1  +  cFitZ); 

end 

line  (radFactCapFit  .*  cFitX,  radFactCapFit  .*  cFitY,  ... 

dim  (2)  *  ones  (size  (cFitX)),  'linestyle',  'color',  back); 

clear  radFactCapFit 
end 

if  showPointGrid 

line  (radFact  (gSqrMaxRow,  gSqrMaxCol)  .*  dirCosXMtx  (gSqrMaxRow,  gSqrMaxCol),  ... 

radFact  (gSqrMaxRow,  gSqrMaxCol)  . *  dirCosYMtx  (gSqrMaxRow,  gSqrMaxCol),  ... 
dim  (2),  'linestyle',  'none',  'marker',  'color',  back); 

end 

if  showPointUnc 
switch  proj 
case  1 

pointRadFact  «  1  /  sqrt  (1  +  pz)  ; 

Jacob  -  (l+pz+px^/pz  px*py/pz  V.  Jacobian  of 

py* px/pz  l  +  pz+pyA2/pz)  /  (2*  (1+pz)  *  (3/2) ) ;  "■  projection 

case  2 

pointRadFact  =  1  /  (1  +  pz) ; 

Jacob  *=  ( l+pz+pxA2/pz  px*py/pz  Jacobian  of 

py*px/pz  l+pz+pyA2/pz]  /  ( 1+pz)  A2 ;  ‘l  projection 

end 

pVarProj  -  Jacob  *  pVar  *  Jacob*;  V.  covariance  matrix  in  this  projection 
IpEigVec,  pEigVall  ...  V.  diagonalize:  pVar 

-  eig  (pVarProj);  *'■  -  pEigVec  *  pEigVal  *  pEigVec* 

pPrAx  -  2  *  pEigVec  ...  V.  scale  principal  axes  (columns  of 

*  sqrt  (pEigVal);  pEigVec)  to  2  standard  deviations 
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line  (pointRadFact  *  px  +  1-1  1]'  *  pPrAx  (1,  :),  ... 
pointRadFact  *  py  +  [-1  11'  *  pPrAx  (2,  :),  — 
dim  (2)  *  111],  ... 

'linestyle',  *-',  'color',  back); 
clear  pointRadFact  Jacob  pVarProj  pEigVec  pEigVal  pPrAx 
elseif  showPoint 
switch  proj 
case  1 

pointRadFact  *  1  /  sqrt  (1  +  pz) ; 
case  2 

pointRadFact  =  1  /  (1  +  pz) ; 

end 

line  (pointRadFact  *  px,  pointRadFact  *  py,  dim  (2),  ... 

'linestyle*,  'none',  'marker',  'color',  back) 

clear  pointRadFact 
end 

if  showSidelobeGrid 

line  (radFact  (sllndx)  .*  dirCosXMtx  (sllndx),  ... 
radFact  (sllndx)  .*  dirCosYMtx  (sllndx),  ... 
dim  (2)  *  ones  (size  (sllndx)),  ... 

'linestyle*,  'none*,  'marker*,  's',  'color',  back); 

end 

if  showSidelobeNrst 
switch  proj 
case  1 

radFactSlNrst  -  1  /  sqrt  (1  +  slNrstVec  (3)); 
case  2 

radFactSlNrst  -  1  /  (1  +  slNrstVec  (3)); 


end 

line  (radFactSlNrst  .+  slNrstVec  (1),  ... 
radFactSlNrst  . *  slNrstVec  (2),  ... 

dim  (2),  'linestyle',  'none',  'marker',  ’  +  ',  'color',  back); 
clear  radFactSlNrst 
end 

if  showSidelobeLgst 
switch  proj 
case  1 

radFactSlLgst  «=  1  /  sqrt  (1  +  slLgstVec  (3)); 
case  2 

radFactSlLgst  -  1  /  (1  +  slLgstVec  (3)); 

end 

line  (radFactSlLgst  .*  slLgstVec  (1),  ... 
radFactSlLgst  .*  slLgstVec  (2),  ... 

dim  (2),  'linestyle',  'none',  'marker',  'x',  'color',  back); 
clear  radFactSlLgst 
end 

end  ‘4  annotations 
*4 

‘4  Arrange  graphics  properties 


axesPat  -  gca; 

rotate  (get  (axesPat,  'children'),  ...  4  rotate  plotted  objects  to 

(0  0  1),  azmthOf f st  /  rpd) ;  %  compensate  for  azmthOffst 

set  (hSurf,  'edgecolor*,  'none'); 
set  (hSurf,  'facecolor',  faceColor) ; 
colormap  (cmap) ; 
set  (axesPat,  'dim',  dim); 

if  cbarVert  V,  position  the  colorbar 

axesCbar  =  colorbar  ('vert'); 

set  (axesCbar,  - 

•units',  'normalized*,  ... 

'position',  ( (l+0.2*cbarSize)/(l+cbarSize)  0.05  0.4*cbarSize/ (1+cbarSize)  0.9]); 
set  (axesPat,  ... 

'units',  'normalized', 

'position' ,  [0  0  1/ (1+cbarSize)  1]); 

else 

axesCbar  =  colorbar  ('horiz'); 

set  (axesCbar,  - 

'units',  'normalized',  - 

♦position*,  (0.05  0 . 3*cbarSize/ (1+cbarSize)  0.9  0 . 5*cbarSize/ (1+cbarSize) ]) ; 
set  (axesPat,  ... 

'units',  'normalized',  ... 

'position',  (0  cbarSize/ (1+cbarSize)  1  1/ (1+cbarSize) ]) ; 

end 

set  (figPat,  'children',  V.  put  color  bar  in  front  of  pattern  but  let 

[axesCbar  axesPat]');  V.  axesPat  remain  the  current  axis 


set  (axesPat,  . . . 

’ xlim’ ,  xylim,  . . . 

'ylim' ,  xylim,  . . . 

' zlim' ,  zlim  ,  ... 

'dataAspectRatio' ,  diff  ([xylim’  xylim'  zlim']),  ... 
'visible',  'off',  ... 

•view',  [viewAz/rpd  viewEl/rpd] ) ; 
if  proj  ==  3 

cameraDist  =  norm  (  ... 

get  (axesPat,  '  cameraPosition' )  ...  V. 

-  get  (axesPat,  ' cameraTarget ' ) ,  2) ; 
set  (axesPat,  *  caraeraViewAngle* ,  ...  '* 

2  *  atan  (float  *  diff  (zlim)  ...  4 

/  (cameraDist  *  diff  (xylim)))  /  rpd);  V, 


distance  from  camera 
to  the  surface 

set  view  angle  to 
exactly  span  the 
unit  sphere 


end 

if  (pointZoom  0)  &  peakVisb 
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pxRot  -  cos  (azmthOffst)  *  px  -  sin  (azmthOffst)  * 
pyRot  -  sin  (azmthOffst)  *  px  +  cos  (azmthOffst)  * 
if  proj  *==  1 

pointRadFact  -  1  /  sqrt  (1  +  pz) ; 


set  (axesPat,  ’xlira’. 

pointRadFact 

*  pxRot  + 

(-1 

* ylim' , 

pointRadFact 

*  pyRot  + 

{-1 

elseif  proj  «  2 

pointRadFact  *1/(1 

+  pz)  ; 

set  (axesPat,  ’xlim'. 

pointRadFact 

*  pxRot  + 

[-1 

' ylim' , 

pointRadFact 

*  pyRot  + 

[-1 

py; 

py; 


i]  / 
i]  / 


i]  / 
i]  / 


elseif  proj  3 

set  (axesPat,  ’cameraViewAngle * ,  — 

2  *  atan  (float  /  (cameraDist  *  pointZoom) )  /  rpd) ; 


V.  compensate  for 
4  azmthOffst 


pointZoom, 
pointZoom) ; 


pointZoom,  . . . 
pointZoom) ; 


end 

clear  pxRot  pyRot  pointRadFact 
end  V,  zoom 

clear  dirCosXMtxSurf  dirCosYMtxSurf  dirCosZMtxSurf ; 
clear  fore  back  gridPolar  gridAzrath  gridAzmthSmth; 
clear  radFactSurf  gridFact  viewAz  viewEl  float  radFact; 


end 


drawnow; 

end  4  loop  over  independent  variable 
clear  indx 

4  Print  warnings,  if  any 

V, 

if  any  (numAcc  <  numRlz) 

disp  (’Warning:  at  least  one  realization  could  not  be  fully*); 
disp  (’  analyzed;  inspect  the  matfiles  in  the  current  directory.*); 
end 

clear  cbarSize 

clear  lx  ly  mx  my  nx  ny; 

clear  numLMX  numLMY  numLMNX  numLMNY; 

clear  tlndx  rlzNum; 

V;  end  of  program 

V,  Suggestions  for  improvements 
% 

V.  Implement  non-uniform  excitation  magnitudes: 

4 

4  As  needed,  non-uniform  illumination  weights  may  be  coded 
%  straightforwardly. 

V. 

V.  Implement  element  factor  and  improve  integration  of  radiated  powers: 

% 

V,  Currently  the  element  factor  is  implicitly  coded  as  one  everywhere, 

%  so  that  the  elements  radiate  uniformly  into  the  hemisphere.  A  more 

%  realistic  element  factor  is  cos  polar  (1),  being  one  at  broadside 

4  and  zero  at  grazing.  Note,  however,  that  in  integrating  the  data 

4  over  solid  angles  to  calculate  radiated  powers,  we  divide  the  data 
4  by  cos  polar.  Performing  this  multiplication  and  division  in 

4  succession  could  yield  invalid  data  near  grazing,  where  cos  polar  is 

4  small.  Obviously,  this  difficulty  could  be  avoided  by  maintaining 

4  the  unsealed  data  for  use  in  the  integration  while  using  the  scaled 

4  data  for  all  other  processing.  The  calculations  of  the  pointing 
4  vector  from  the  excitation  phases  and  of  the  peak  power  should  also 
4  account  for  the  element  factor. 

4 

4  More  generally,  one  might  wish  to  apply  an  arbitrary  element  factor. 
4  If  it  is  small  near  grazing,  the  difficulty  described  above 
4  persists.  One  solution  is  to  specify  the  element  factor  relative  to 
%  COS  polar;  the  user  ensures  that  all  values  of  that  ratio  are 
4  reasonable.  The  data  used  in  the  integration  are  scaled  by  the 

4  ratio,  while  the  data  used  elsewhere  are  scaled  by  both  the  ratio 

4  and  cos  polar.  One  could  also  allow  element  factors  defined  over 

4  the  entire  sphere,  including  radiation  into  z  <  0,  perhaps  by 
4  storing  the  back  radiation  pattern  in  a  second  g  matrix  and 
4  modifying  the  analysis  routines.-, 

4  Broadening  our  perspective,  we  note  that  these  problems  are 

4  ultimately  due  to  the  integration  algorithm,  in  that  it  blindly 

4  applies  a  Jacobian  that  diverges  at  grazing.  One  consequence  is 
4  that  data  cells  whose  centers  are  just  inside  the  border  of  visible 

4  space  contribute  their  entire  value  scaled  by  a  large  Jacobian, 

4  whereas  those  whose  centers  are  just  outside  contribute  nothing. 

4  More  correctly,  both  should  contribute  about  half  of  their  value 

4  scaled  according  to  some  average  location  of  the  contributing 
4  region.  Cells  almost  entirely  outside  of  visible  space  should 
4  contribute  very  little.  A  better  algorithm  would  reproduce  this 
4  behavior. 

4  Allow  pattern  cuts: 

4  Often  one  is  interested  in  a  cut  of  the  pattern  along  a  path  on  the 
4  unit  hemisphere,  such  as  cuts  through  the  main  beam  along  azimuth 
4  and  elevation  curves  or  along  the  great  circles  of  maximum  and 
4  minimum  beam  width.  If  only  a  graphical  presentation  is  desired, 

4  the  cut  could  be  interpolated  from  the  pattern.  If  an  analysis  is 
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also  desired,  specialized  routines  would  probably  be  required,  as 
\\  the  current  routines  expect  data  over  two  directional  coordinates . 

V.  Consider  alternate  calculation  of  beam  widths  and  roll: 


■*  The  beam  widths  and  roll  are  derived  from  an  ellipse  fitted  to  a 
%  level  contour  of  the  main  beam.  Currently  the  fit  is  performed  in 
the  boresight-centered  stereographic  projection  regardless  of  the 
1  pointing  vector.  This  projection  preserves  the  orientation  and 
\\  orthogonality  of  the  major  and  minor  axes  of  the  ellipse.  However, 

•;  because  scale  in  the  stereographic  projection  increases  away  from 

boresight,  of f-boresight  contours  are  expanded  toward  the  edge  of 
V.  the  projection.  Although  scale  is  uniform  in  all  directions  for  an 
infinitesimal  region,  radial  scale  is  exaggerated  relative  to 
V,  azimuthal  scale  for  a  finite  region.  For  a  broad  beam  off 

S  boresight,  the  distortion  of  beam  width  may  be  significant. 


An  improvement  might  be  achieved  with  two  changes.  First,  center 
the  projection  on  the  pointing  vector,  assuming  that  the  contour 
will  be  found  to  lie  centered  on  the  pointing  vector  as  well. 

Second,  instead  of  the  stereographic  projection,  use  the  azimuthal 
equidistant  projection,  for  which  distances  measured  on  a  line 
passing  through  the  center  of  the  projection  are  true.  Regardless 
of  beam  width,  the  lengths  of  the  major  and  minor  axes  of  the 
projected  ellipse  will  equal  those  of  the  ellipse  on  the  hemisphere, 
provided  that  the  center  of  the  ellipse  is  also  the  center  of  the 
projection. 

In  centering  the  projection  on  the  pointing  vector,  one  is  free  to 
choose  the  orientation  of  the  projection  relative  to  the  spherical 
coordinate  system,  and  a  judicious  choice  of  this  angle  facilitates 
calculation  of  the  beam's  roll  angle.  Construct  at  boresight  on  the 
hemisphere  two  tangent  axes  u  and  v;  let  positive  u  be  directed 
toward  spherical  azimuth  zero  and  positive  v  toward  pi/2.  Next, 
construct  the  arc  connecting  boresight  with  the  pointing  vector. 
Translate  the  u-v  origin  and  system  along  this  arc  without  rotation 
in  the  plane  locally  tangent  to  the  hemisphere,  so  that  u,  v,  and 
the  arc  maintain  the  same  local  orientation.  If  the  ellipse  is 
centered  on  the  pointing  vector,  the  roll  angle  will  be  the  angle 
from  the  positive  u  axis  to  the  major  axis. 

Allow  linear  arrays; 


Vi 

V.  The  current  analysis  routines  cannot  handle  linear  arrays  because  of 
%  several  incompatibilities.  For  example,  the  main  beam  of 
%  one-dimensional  arrays  with  isotropic  element  patterns  is  a  cone 
about  the  axis  of  the  array,  so  the  direction  of  radiation  is 
%  specified  by  a  single  number  —  the  angle  between  the  axis  and  the 
%  cone.  However,  the  current  analysis  routines  seek  a  two-parameter 
%  specification  of  the  main  beam  direction  and  will  generally  fail. 

%  Likewise,  the  beam  width  is  described  by  one  number,  but  the  current 

%  routines  seek  two  parameters.  Furthermore,  analyses  that  depend  on 
%  these  values  (such  as  determining  the  proximity  of  sidelobes)  will 
4  also  fail.  If  linear  arrays  are  if  interest,  the  analysis  routines 
%  must  be  expanded.  One  might  also  add  graphics  routines  tailored  to 
'{,  one-dimensional  arrays. 


V.  Note  that  a  non-isotropic  element  pattern  will  generally  produce  a 
'4.  variation  in  the  direction  orthogonal  to  the  main  beam  contour, 

‘4  allowing  the  analysis  routines  to  proceed.  Results  thereby  obtained 

%  should  be  interpreted  accordingly.  (Using  a  non-isotropic  element 

%  pattern  will  not  benefit  the  routine  that  locates  the  pointing 
%  vector  from  the  excitation  phases.  One  might  ignore  its  results 
when  the  final  pointing  vector  is  calculated,  using  only  the 
l  pointing  vector  obtained  from  the  Fourier  transform.) 

?. 

\*  [1]  R.  Tang  and  R.  W.  Burns,  "Phased  Arrays,"  in  _Antenna 

V.  Engineering  Handbook^,  3rd  ed.,  R.  C.  Johnson,  Ed.  New  York,  NY: 

V,  McGraw-Hill,  1993,  Ch.  20,  Sec.  3. 
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